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Abstract 

A method for calculating unsteady flows in cascades is presented. The model, which 
is based on the linearized unsteady Euler equations, accounts for blade loading, shock 
motion, wake motion, and blade geometry. The mean flow through the cascade is 
determined by solving the full nonlinear Euler equations. Assuming the unsteadiness in 
the flow is small, then the Euler equations are linearized about the mean flow to obtain a 
set of linear variable coefficient equations which describe the small amplitude, harmonic 
motion of the flow. These equations are discretized on a computational grid via a 
finite volume operator and solved directly subject to an appropriate set of linearized 
boundary conditions. The steady flow, which must be calculated prior to the unsteady 
flow, is found via a Newton iteration procedure. Using this procedure, the nonlinear 
steady Euler equations are solved as a series of linear problems, each very similar to 
the unsteady linear problem. The similarity of the steady and unsteady solvers greatly 
reduces the effort which is needed to develop the two codes. 

An important feature of the analysis is the use of shock fitting to model steady 
and unsteady shocks. Use of the Euler equations with the unsteady Rankine-Hugoniot 
shock jump conditions correctly models the generation of steady and unsteady entropy 
and vorticity at shocks. In particular, the low frequency shock displacement is correctly 
predicted. 

Results of this method are presented for a variety of test cases. Predicted unsteady 
transonic flows in channels are compared to full nonlinear Euler solutions obtained 
using time-accurate, time-marching methods. The agreement between the two meth- 
ods is excellent for small to moderate levels of flow unsteadiness. The method is also 
used to predict unsteady flows in cascades due to blade motion (the flutter problem) 
and incoming disturbances (the gust response problem). Comparison of the predicted 
unsteady flow to other semi-analytical and numerical methods gives good agreement. 
The linearized Euler method requires substantially less computational effort than the 
time-marching procedures, making the present method useful for aeroelastic analyses. 


2 


Acknowledgments 


The primary funding for this research was provided by a grant from NASA Lewis Re- 
search Center (Grant Number NSG 3079) which was monitored by Dr. John Adam- 
czyk, Mr. Donald Boldman, and Dr. Dan Hoyniak. Additional support was provided^ 
by the Aircraft Engine Business Group of the General Electric Company, the Alison 
Gas Turbine Operations of the General Motors Corporation, and the Lockheed-Georgia 
Company. The author gratefully acknowledges the generous support of the Fannie and 
John Hertz Foundation. 


3 


Contents 


1 Introduction 6 

2 Analysis of Quasi-One-Dimensional Flows 14 

2.1 Linearized Quasi-One-Dimensional Euler Equations 15 

2.2 Solution of the Steady E..*ier Equations 18 

2.2.1 Newton Iteration Equations 18 

2.2.2 Discretization of the Newton Iteration Equations 19 

2.2.3 Boundary Conditions 20 

2.2.4 Assembly of the Discrete Equations 21 

2.2.5 Subsonic Flow in a Converging-Diverging Channel 26 

2.2.6 Steady Shock Fitting 26 

2.2.7 Steady Transonic Channel Flow 32 

2.3 Solution of the Linearized Unsteady Euler Equations 36 

2.3.1 Discretization of the Linearized Euler Equations 39 

2.3.2 Unsteady Channel Flow Boundary Conditions 39 

2.3.3 Assembly of the Discrete Equations 40 

2.3.4 Unsteady Subsonic Results 40 

2.3.5 Unsteady Shock Fitting 42 

2.3.6 Assembly of Unsteady Transonic Equations 45 

2.3.7 Unsteady Transonic Results 45 

2.4 Summary 48 

3 Advantages and Limitations of the Linearized Euler Analysis 49 

3.1 Nonlinearities in Unsteady Flow 50 

3.1.1 Fourier Analysis of Unsteady Flows 51 

3.1.2 Numerical Examples 55 

3.2 Comparison to Full Potential Theory 65 

3.3 Summary 75 

4 Fundamentals of Two-Dimensional Linearized Euler Theory 76 

4.1 The Conservation Laws 77 

4.1.1 Integral Form of the Conservation Laws 77 

4.1.2 Integral Form of the Conservation Laws for Deforming Volumes . 79 

4.1.3 Differential Form of the Conservation Laws 80 

4.2 Linearized Euler Equations 82 


4 


4.2.1 Field Equation Linearization 83 

4.3 Solution of the Steady Euler Equations 86 

4.3.1 The Linearized Newton Iteration Equations 86 

4.3.2 Discretization of the Newton Iteration Equations 87 

4.3.3 Steady Flow Boundary Conditions 92 

4.3.4 Matrix Equation Structure 94 

4.4 Unsteady Flow Calculations 97 

4.4.1 Linearized Unsteady Euler Equation Discretization 98 

4.4.2 Boundary Conditions of the Unsteady Flow Solution 99 

4.4.3 Matrix Equations 103 

4.5 Summary 103 

5 Extensions of Linearized Euler Theory to Cascade Flow 104 

5.1 Moving Airfoil Boundaries 106 

5.2 Flow Discontinuities 108 

5.2.1 Unsteady Shock Jump Conditions 109 

5.2.2 Unsteady Wake Jump Conditions 112 

5.2.3 Newton Iteration Solution Procedure for Flows with Shocks and 

Wakes 114 

5.3 Periodic and Far Field Boundary Conditions 116 

5.4 Surface Pressures 122 

5.4.1 Moving Airfoils 122 

5.4.2 The Shock Impulse 123 

5.5 Summary 123 

6 Results 125 

6.1 The Flow in a Hyperbolic Channel 125 

6.1.1 Accuracy of the Steady Solver 127 

6.1.2 Convergence of Subsonic Solutions 132 

6.1.3 Steady Choked Flow in a Hyperbolic Channel 134 

6.1.4 Effect of Grid Resolution on Shock Location 137 

6.1.5 Unsteady Subsonic Flow in a Hyperbolic Channel 137 

6.1.6 Unsteady Transonic Flow in a Hyperbolic Channel 139 

6.2 Two-Dimensional Transonic Diffuser Flow 146 

6.2.1 Steady Diffuser Flow 147 

6.2.2 Unsteady Transonic Diffuser Flow 149 

6.3 The Gostelow Cascade . 156 

6.3.1 Steady Flow Through the Gostelow Cascade 156 

6.3.2 Unsteady Flow in the Gostelow Cascade 159 

7 Conclusions and Recommendations 177 


5 


Chapter 1 


Introduction 


The flutter and forced response of turbomachinery blading has been and continues to 
be a recurring problem in the development of turbomachinery. Nevertheless, the ability 
to predict these two aeroelastic phenomena has remained illusive. Although the struc- 
tural dynamic modelling of rotors has steadily progressed over the years, the modelling 
of the unsteady aerodynamics of cascades still needs much improvement. The objective 
of this research is to model unsteady transonic flows in cascades using a computational 
scheme based on the linearized Euler equations. To date, much theoretical work has 
been performed in the field of unsteady flow in cascades. Most of it, however, has been 
based on the assumption of isentropic irrotational flow even through shocks. These 
methods do not allow for the production of steady and unsteady entropy and vortic- 
ity across shocks. This weakness in potential theory can, under certain circumstances, 
produce serious errors in the predicted shock motion and unsteady blade loading. A 
new method based on the linearized Euler equations will be presented in this report. 
The method correctly accounts for the production of steady and unsteady entropy and 
vorticity across shocks and therefore predicts unsteady transonic flows more accurately. 

Theoretical research into unsteady flows in cascades has steadily progressed over 
the past three decades. The methods used can be roughly divided into three types: 
Analytical methods, semi-analytical methods, and numerical methods. In analytical 
methods, the partial differential equations which govern the unsteady flow in a cascade 
are solved exactly using applied mathematical techniques. Although closed-form solu- 
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tions can sometimes be obtained, the resulting complicated expressions must usually 
be evaluated numerically to understand them. Furthermore, these solutions have only 
been found for extremely simplified geometries. The semi-analytical methods use singu- 
larities such as vortices distributed along the surface of the airfoils and wakes to model 
the steady and unsteady flow. In fact, the distinction between the analytical and semi- 
analytical methods is small. In practice, both usually reduce the problem to integral 
equations involving kernel functions. Because of the complexity in determining these 
kernel functions for complicated flows, the mean flow through the cascade is usually 
assumed to be uniform or nearly uniform. Hence, the flow model often bears little re- 
semblance to flows found in actual turbomachinery. The third approach is to discretize 
the field equations which describe the unsteady flow and solve them numerically. This 
computational approach has the advantage that arbitrary airfoil shapes can be ana- 
lyzed. Furthermore, flow features such as blade thickness, blade loading, and moving 
shocks and wakes are accounted for. This ability to analyze fairly general geometries 
is obtained only with considerable computational effort. The computational methods 
can be subdivided further into two categories: time-marching and linearized harmonic. 
With time marching methods, the unsteady flow field is found as a function of time 
by marching a flow simulation in time. This approach, although straight forward, is 
computationally expensive as many time steps must be taken to reach a periodic state. 
In the linearized harmonic methods, the nonlinear field equations are linearized about 
some nominal mean flow, resulting in a set of linear variable coefficient field equations 
for the unsteady flow perturbations. If the additional simplifying assumption is made 
that the unsteady flow is harmonic in time, these equations are reduced to a set of time- 
independent partial differential equations for the complex amplitudes of the unsteady 
flow properties. These linearized harmonic field equations are discretized and solved 
directly. This is the approach taken in the present research. 

Previous Work on Unsteady Aerodynamics 

Whitehead [l] studied the unsteady flow of an incompressible fluid through a cascade 
of vibrating flat plate airfoils. In this model it is assumed that the mean flow is uniform 
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and undeflected by the cascade. Unsteady vorticity is distributed along the airfoils and 
shed vorticity is convected along the wake. The bound vorticity is distributed such 
that the velocity due to the vorticity is equal to the upwash velocity on the airfoil 
surfaces and the Kutta condition is satisfied. This model, however, failed to predict 
the experimentally observered bending flutter. Later Whitehead [2] extended his model 
to that of a cascade of flat plat airfoils for which the mean flow was not uniform, but 
was deflected or turned by the cascade. This new model did predict bending flutter 
and showed the importance of steady loading on the unsteady pressure distributions 
on vibrating airfoils. Atassi and Akai [3,4] used a model in which point singularities 
were distributed along the surface of a two-dimensional airfoil surface of finite camber 
and thickness. With this model they studied steady and unsteady incompressible flows 
through a cascade of thick cambered airfoils. Their results, like Whitehead’s, showed 
the importance of the steady blade loading on the unsteady flow in the cascade. 

Lane and Friedman [5], Whitehead [6], and Smith [7] all analyzed flat plate cas- 
cades vibrating in a uniform subsonic flow. With compressibility there is the added 
complication of acoustic modes and acoustic resonance. Under certain circumstances, 
acoustic waves will propagate away from the rotor unattenuated, while in other cases, 
the waves will be cut-off. Smith [7] performed experiments to validate his model. For 
unloaded cascades, the theory is successful in predicting the cut-off behavior as well as 
the amplitude of the acoustic waves. For steadily loaded cascades, the amplitudes of 
the acoustic waves were not predicted as well. 

Several investigators [8,9,10,11,12] have studied the problem of a cascade of vibrating 
flat plate airfoils in a uniform supersonic flow which is axially subsonic. Although 
the analysis techniques are differ, all these investigators studied the same governing 
differential equations and boundary conditions. Hence, it is not surprising that all 
the results from these models are essentially the same. These models indicate that 
bending flutter will not occur at the reduced frequencies and Mach numbers at which real 
compressors actually exhibit bending flutter. Bendiksen [13] later used a perturbation 
scheme to calculate the effects due to small amounts of thickness, camber, and angle of 
attack as well as shock motion, and demonstrated the important role of shock motion 
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in flutter prediction. Goldstein, Braun, and Adamczyk [14] have proposed a model in 
which the steady flow is everywhere parallel but has a strong nonisentropic in-passage 
shock. Hence, although there is no turning, the blade row diffuses the flow due to the 
presence of a shock. That is to say that the blades carry a steady load due to the pressure 
rise across the shock. Further, in this theory the unsteady shock motion is correctly 
accounted for as well as the vorticity and entropy generation at the shock. Using this 
model, bending flutter is predicted for reduced frequencies (based on semichord) below 
about 0.3. This result is more in keeping with experimentally observed bending flutter 
in compressors operating at high back pressures. 

The lesson learned from these analytical and semi-analytical models is clear. To 
model accurately the flow in a cascade, one must include the effects due to the steady 
blade loading. This means that the camber and thickness of the airfoil need to be 
accounted for, as well as the presence of shocks in the flow. Except for the incompressible 
case, this forces the investigator to turn to numerical methods because it is not possible 
to analytically solve linear field equations with complicated variable coefficients. To 
investigate subsonic flows, Verdon and Caspar [15] numerically determined the mean 
flow about a cascade of airfoils in subsonic flow using the full potential equations. The 
full potential equations were then linearized about the mean flow and solved numerically 
using a finite difference approximation. The computed unsteady pressure distribution 
agreed well with tests conducted by Carta [16] on a linear cascade of oscillating airfoils. 
Whitehead and Grant [17] performed a similar analysis but used finite elements to 
discretize the field equations. Later, Verdon and Caspar [18] extended their model to 
handle transonic flows. To accurately model the effect of shock motion, shock fitting 
was used. Whitehead [19] also extented his model to the transonic regime but used 
shock capturing to model the shock motion. Both models require some form of artificial 
viscosity to stabilize the calculations in the supersonic regions. 

Several investigators have used time-marching schemes to analyze unsteady flows 
in turbomachinery. Ni and Sisto [20] used a time-marching Euler method to calculate 
the pressure loads on vibrating flat plate airfoils in compressible flow. Actually, the 
method can be considered a hybrid of the time-marching and harmonic techniques 
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as the quantities which are time-marched are the harmonic flow variables. Recently, 
Giles [21,22] has developed an Euler time-marching method based on Ni’s Lax-Wendroff 
method [23] to study wake/rotor interaction. The ability to handle unequal rotor/stator 
pitches is accomplished by inclining the computational plane in time. This allows the 
calculations to be performed in a single blade passage. The code has the advantage 
that nonlinear, anharmonic flows can be modelled. However, for numerical stability, the 
time integration step size used in the integration must be fairly small. The mavim nm 
allowable time step is governed by the so-called CFL number. Because of this CFL 
number restriction, the size of the time steps are limited roughly by the size of the 
smallest cell in the domain. Hence, for accurate resolution of the flow, the required 
computational times become quite large. 

Description of the Linearized Euler Analysis 

In reality, both steady and unsteady flows in turbomachinery are extremely complicated. 
The fluid is viscous and heat-conducting and is most accurately described by the Navier- 
Stokes equations. However, if the Reynolds number is sufficiently high and the Prandtl 
number is order unity, and separation doesn’t occur, the viscous and heat transfer effects 
are confined to narrow regions near the airfoil surfaces and in the wakes. Under these 
circumstances, the Euler equations are a good approximation to the behavior of the flow. 
The unsteady Euler equations are the starting point of the linearized Euler analysis. 
The Euler equations correctly account for the production of entropy and vorticity at 
shocks and therefore are the natural choice for calculating unsteady transonic flows. 

The method is divided into two mains parts. The first step is to determine the steady 
flow. In order to fully account for such effects as blade thickness, loading, and strong 
in-passage shocks, the unsteady Euler equations should be linearized about the correct 
mean flow. So despite the fact that we are primarily interested in the unsteady behavior 
of the fluid, the steady flow will also have to be calculated. Because of the periodicity 
in a cascade, the flow cm be computed in a single blade passage. Within this passage, 
the steady Euler equations are discretized on a computational grid using a conservative 
finite volume operator. The discretized equations are solved subject to an appropriate 
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set of boundary conditions using a Newton iteration procedure. Using this procedure, 
the nonlinear Euler equations are transformed into a series of linear equations [24,25]. 
These are solved directly using Gaussian elimination. The Newton iteration procedure 
is extremely efficient and usually converges within five to ten iterations. 

Once the steady solution has been obtained, the unsteady flow can be calculated. 
The nonlinear time-dependent Euler equations are linearized about the steady solution 
to obtain the linearized unsteady equations. These equations are linear with variable co- 
efficients and describe the small disturbance behavior of the flow. Since many unsteady 
flows of interest are periodic in time, the unsteady flow is assumed to be harmonic 
in time. Under this assumption, the explicit time dependency is eliminated from the 
problem. As with the steady Euler equations, the unsteady flow is computed in a single 
blade passage. The equations are discretized on a computational grid using a finite 
volume operator. 

The application of the unsteady boundary conditions is a crucial step in the for- 
mulation of the problem. For example, because the computational domain is finite, an 
artificial boundary must be placed in front of the cascade. Waves which propagate from 
inside the blade passage must pass through this boundary without being reflected if the 
unsteady flow is to be correctly predicted. So-called nonreflecting boundary conditions 
are applied by use of a characteristic wave analysis. 

An important feature of the present method is that steady and unsteady shocks 
and wakes are fitted rather than captured. Most CFD codes use shock capturing to 
model shocks. This usually means that some artificial viscosity is added to the scheme. 
The resulting captured shock is not a sharp discontinuity but is smeared over several 
grid points. The advantage to this approach is that shocks are automatically captured 
where ever they happen to be, but increased grid resolution is needed in the area of the 
shock to produce a good approximation to a jump discontinuity. Using shock fitting, 
the position of the shock is modelled explicitly with unsteady shock -jump conditions 
applied across the shock [26,27] . These conditions, like the Euler equations themselves, 
can be linearized to obtain a set of linearized shock jump conditions. Furthermore, 
the jump conditions fully account for the production of unsteady entropy and vorticity 
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through shocks. This is a major improvement over the current linearized full potential 
theories. In particular, at low reduced frequencies, the shock displacement predicted by 
the linearized potential theory is too large. 

Outline 

In Chapter 2 of this report, the basic method for solving quasi-one-dimensional steady 
and unsteady flows is presented. The purpose is to introduce the fundamental ideas 
without some of the complicating details required to analyze two-dimensional flows. 
Furthermore, the quasi-one-dimensional model problems presented in this chapter are 
ideal test vehicles for validating the basic methodology. In this chapter, all of the funda- 
mental steps needed to analyze steady and unsteady flows are presented. They include 
the linearization of the nonlinear Euler equations, the discretization of the linearized 
Euler equations, the linearization and discretization of unsteady nonlinear boundary 
conditions, and the implementation of the direct solver used to solve these equations. 
Also discussed is the method of steady and unsteady shock fitting in one-dimension. 
Finally, the linearized Euler method is used to analyze model unsteady flow problems. 
These results are compared to those from a time-marching Euler code to demonstrate 
the validity of the linearized Euler method. 

In Chapter 3, some of the advantages and limitations of the linearized Euler method 
are discussed. Using an existing time marching algorithm to solve quasi-one-dimensional 
model problems, the effects of nonlinearities are investigated. As will be shown, the 
nonlinear effects enter mainly through higher harmonics. The fundamental harmonic 
component is accurately represented by the linearized theory for small to moderate 
levels of unsteadiness. Also discussed are the advantages the linearized Euler method 
offers over linearized full potential methods. The latter represent the current state 
of the art in linearized theories for solving for unsteady flows through cascades. The 
fundamental assumption made in the potential methods is that the flow is isentropic and 
irrotational. Hence, although the methods are efficient, they do not correctly account 
for the unsteady shock behavior in transonic flows. This is illustrated by examining 
the linearized Euler and linearized potential solutions for simple quasi-one-dimensional 
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model problems. 

In Chapter 4, the theory and numerical method for solving two-dimensional chan- 
nel flow problems is presented. As in the one-dimensional problem, this process can 
be broken down into the linearization of the Euler equations, the linearization of the 
unsteady boundary conditions, the discretization of the field and boundary-condition 
equations, and the numerical solution of the discretized equations. In Chapter 5, the 
discussion of the two-dimensional theory and numerical method continues. Discussed 
are the extensions of the basic theory needed to calculate cascade flows. These include 
discussions of steady and unsteady shock fitting, wake fitting, moving airfoil boundary 
conditions, and nonreflecting far-field boundary conditions. 

In Chapter 6, a variety of two-dimensional steady and unsteady flow calculations 
are presented. These include steady subsonic and transonic channel flows, unsteady 
subsonic and transonic channel flows, and steady and unsteady cascade flows. The 
unsteady cascade flows presented are of two types. The flutter problem, where the 
motion of the blades are prescribed and the unsteady blade loads are to be calculated, 
and the gust response problem, where an incoming gust is specified and the unsteady 
airloads on the stationary blades are to be calculated. 

Finally in Chapter 7, some conclusions are drawn about the utility of the present 
method, and suggestions are given for areas of future research needed to improve the 
present method. 
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Chapter 2 


Analysis of 

Quasi- One-Dimensional Flows 


In this chapter, the linearized Euler analysis is applied to quasi-one-dimensional 
channel flow problems. By studying these relatively simple model problems, the funda- 
mental linearized Euler analysis can be illustrated without some of the complications 
which arise in two-dimensional problems. In Section 2.1 the linearization of the Euler 
equations is examined. This linearization reduces the nonlinear time-dependent Euler 
equations to a set of linear variable-coefficient partial differential equations. The lin- 
ear equations describe the small disturbance behavior of the unsteady flow. Because 
the Euler equations are linearized about a steady flow solution, this solution must be 
determined before the analysis of the linearized unsteady flow can be performed. In 
Section 2.2 the method for computing steady Euler flows is presented. The nonlinear 
Euler equations are solved using a Newton iteration technique. This method is not only 
computationally efficient, but also reduces the nonlinear problem to a series of linear 
problems each similar to the linearized unsteady Euler equations. This greatly reduces 
the development effort required to write both a steady and unsteady code since many 
the components will be nearly identical. Some sample subsonic and transonic quasi- 
one-dimensional flows will be calculated using this approach. In Section 2.3 the method 
for solving the linearized Euler equations will be discussed. As with the steady prob- 
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lem, the steps to be taken are the discretization of the field equations, and boundary 
conditions, the assembly of the equations into matrix form, and the solution of these 
linear equations. 

2.1 Linearized Quasi-One-Dimensional Euler Equations 

Consider the unsteady Euler equations for flow in a quasi-one-dimensional channel 
with area variation A(x) as shown in Figure 2.1. The mass, momentum, and energy 
conservation equations are 

3. AV+ ‘L AF + A JL T=0 ( 2 .,) 

where 



where p , u, p, e<)» and h o are the static density, velocity, static pressure, toted interned 
energy per unit volume, and toted specific enthalpy, respectively. By definition, the total 
specific enthedpy is given by 

ho = — (2-2) 

P 

The toted interned energy is the sum of the intrinsic energy plus the kinetic energy. For 
an ideed gas with constant specific heats, the intrinsic energy is c v T = J. Under 
this assumption, Equation 2.2 can be rewritten as 

< 2 - s > 

Elimin ation of the total enthalpy and energy in favor of the three other primitive vari- 
ables gives 

p pu 0 

U = pu F = pu* P = p 

. ^ + \f> ul . . + \ P u * . . 0 . 

This form of the Euler equations will be used throughout, although the methods used 
here could be extended to consider a non-ideal gas if desired. These unsteady Euler 


15 


equations are strictly valid for inviscid adiabatic flow of a perfect gas. If the Reynolds 
number is sufficiently large and no separation occurs, however, they provide an adequate 
representation of the flow outside thin viscous boundary layers and wakes. 

Note in particular that the unsteady Euler equations are nonlinear in the primi- 
tive variables p, u, and p. An exact solution would have to take into account these 
nonlinearities. One approach is to simulate the flow using a time marching simulation. 
Given enough spatial and temporal resolution and computer time, one could find the 
full nonlinear time-dependent solution. This approach, while straightforward, would 
be computationally expensive, and for many problems of interest, the effects of the 
nonlinearities are small and can be ignored. 

The approach taken in the present research is to linearize the Euler equations about 
some nominal mean flow. If the unsteady part of the flow is small compared to the 
mean flow values, then this is a reasonable approximation to the full nonlinear Euler 
equations. The linearized equations could also be solved by time-marching. However, 
many unsteady flows of interest are periodic in time, especially in the area of aeroelas- 
ticity. Hence, the further simplification can be introduced that the unsteady part of the 
flow is harmonic in time. This assumption removes the explicit time dependency from 
the governing equations. Hence, the unsteady problem is reduced to solving a set of 
linear time-independent ordinary differential equations. 

To start the linearization process, the flow can be expressed as the sum of two parts: 
a mean or steady flow component which satisfies the steady nonlinear Euler equations, 
and a small unsteady perturbation which is harmonic in time, i.e., 


p{x, y, t) = p{x, y) + p(x, y)e’ ut 

(2.4) 

u(*,y,t) = U{x, y) + u(x, y)e* ut 

(2.5) 

p(x, y, t) = P{x, y) -1- p(x, y)e jut 

(2.6) 


where the perturbation variables p , u, and p are small compared to the mean-flow coun- 
terparts p, U, and P, and are complex due to the harmonic motion assumption. Note 
the change in notation made to distinguish the different parts of the flow. The steady- 
flow variables are represented by upper case characters or by lower case characters with 
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an overbar, the perturbation variables are represented by lower case characters, and the 
total flow which is the sum of the two is now represented by lower case characters with 
a circumflex. 

Substitution of the foregoing expansion into the one-dimensional Euler equations 
and collection of the zeroth- and first-order terms results in the one-dimensional steady 
Euler and linearized unsteady Euler equations. The zeroth order equations are just the 
original nonlinear Euler equations (Equation 2.1) with the time derivative term set to 
zero. For small unsteady perturbations the mean flow is independent of the unsteady 
part of the flow. The first-order equations are given by 

d d 

jwABiU -I- — ABjU + A— B311 = 0 (2-7) 

where u is the vector of perturbation variables p, u, and p, and Bi, B2, and Bj are 
variable coefficient matrices which depend on the mean flow variables p t U , and P. The 
coefficient matrices are given by 

10 0 
Up 0 

V* A. 

p 0 

2 pU 0 

P+IPU 2 

0 0 0 
0 0 1 
0 0 0 

These unsteady perturbation equations, and their two-dimensional equivalents (to be 
derived in Chapter 4), are the foundation of the analysis used in this work to analyze 
unsteady flows. They are valid so long as the unsteady perturbation variables are 
small compared to the mean-flow variables. In the next chapter, it will be shown using 
numerical experiments that the perturbations equations are a useful approximation to 
the full nonlinear Euler equations for fairly moderate levels of unsteadiness. 



B,= 


dU 

du 


02 du 


U 

U* 
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The general approach for solving unsteady flows is as follows. First, the mean flow 
is calculated using the steady Euler equations. This solution is then used to calculate 
the variable coefficients in the linearized unsteady equations. The linearized Euler 
equations are discretized and solved numerically to obtain the unsteady flow solution. 
The remainder of this chapter is devoted to explaining the details of this approach as 
applied to subsonic and transonic channel flows, as well as to presenting solutions to 
some model problems. 

2.2 Solution of the Steady Euler Equations 

2.2.1 Newton Iteration Equations 

Before the unsteady flow can be calculated, one must first have solved for the steady 
flow. Although the purpose of this research is to better understand unsteady flows and 
to develop a linearized Euler solver, a steady Euler solver is also needed to calculate 
the steady flow which is used as the input to the unsteady Euler solver. Fortunately, 
the steady nonlinear Euler equations can be solved using a Newton iteration proce- 
dure [24,25] which converts the nonlinear ordinary differential equations into a series of 
linear ordinary differential equations each very similar to the linearized Euler problem. 
Again, the flow is assumed to be composed of two components: a base solution which 
is not assumed to satisfy the steady Euler equations, plus a small perturbation quan- 
tity which when added to the base solution will give an improved estimate of the true 
solution. 

p{x, y, t) = p{x , y) + p( x, y) (2.8) 

«(*, y, t) = U[x, y) + u(*, y) (2.9) 

p(*, y , «) = P(*, y) + p(®, y) (2 . 10) 

where the correction perturbation variables p, u, and p are assumed to be small. Substi- 
tution of Equation 2.8-2.10 into the steady Euler equations, again collecting first-order 
perturbation terms, and neglecting higher order effects gives 

+ a £ BsU = “ (s AF + A ti F ) ( 2 “) 
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The terms on the right-hand-side are the residuals of the full nonlinear steady Eu- 
ler equations evaluated from the current estimate of the flow ( p , U , and P) through 
the channel; those on the left-hand-side represent the first-order effect the correction. 
Equation 2.11 is solved numerically to determine the correction perturbations. These 
corrections are then added to the previous estimate of the solution to obtain an improved 
solution. Then, starting with this updated estimate of the flow, the entire process is 
repeated until it converges, although convergence is not guaranteed. When the method 
works, convergence is usually quite fast with a converged solution obtained in five to 
ten iterations. 

Note the similarity of the Newton iteration perturbation equation (Equation 2.11) 
to the linearized unsteady Euler equations (Equation 2.7). The matrices B 2 and Bs in 
Equation 2.11 are the same as those in Equation 2.7. The Newton iteration technique is 
an extremely efficient one for calculating steady one-dimensional and two-dimensional 
flows. Furthermore, because of the similarities to the unsteady problem, the effort re- 
quired to construct the steady and unsteady codes is greatly reduced. Many components 
of the two codes are nearly identical. 

2.2.2 Discretization of the Newton Iteration Equations 

Consider first the discretization of the linearized steady Euler equations (Equa- 
tion 2.11). These equations are discretized on a one-dimensional computational grid as 
shown in Figure 2.1. The grid has 7 nodes and hence 7—1 conservation cells. At each 
of the nodes, the values of the primitive variables are stored. For each of the 7—1 cells, 
the conservation equations are approximated by a central difference operator operating 
about the cell center. Hence, the Newton iteration equations for the steady flow are 
approximated at the ith cell by 

(AB 2 u), + 1 - (AB 2 u)< , (B 3 u) l+1 - (B s u),. _ 

5^ + M Ax“ 

Axi M A“ (212) 

where 

AXg — 2g+l 
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Figure 2.1: Channel geometry and one-dimensional grid use in numerical calculation of 
steady and unsteady flow 

A i+$ = $( A » + A »+i) 

It can be shown that these difference equations are second order accurate. That is to 
say that the difference between the computed solution and the exact solution will be 
0(Ax 2 ). These difference equations provide 3(f— 1) equations for 31 unknowns. Still 
to be specified are the upstream boundary and downstream boundary conditions. The 
boundary conditions, like the field equations must be linearized for use in the Newton 
iteration solver. This linearization will be discussed in the next section. 


2.2.3 Boundary Conditions 


The boundary conditions, like the Euler equations, will, in general, be nonlinear 
functions of the primitive variables. Hence, to apply the boundary conditions in the 
Newton iteration scheme requires that they be linearized. In particular, for subsonic 
inlet flow, two boundary conditions must be specified. Usually the total pressure and 
density are given at the inlet of the channel. The total pressure and density are related 
to the velocity and the static values of pressure and density by 


TP 7 ~ 1 -2 _ 7firo 

p 2 pro 


(2.13) 
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P 1 Pro 

These two inlet boundary conditions are linearized for the Newton iteration procedure 
just as the field equations were, i.e., by expanding the nonlinear boundary conditions 
in a perturbation series and collecting terms of first order to give 


7 7 P t f «\rr iPrO 

-p ~ ~=»P + (7 - 1 )U u = 
p p 2 PTO 


(2.15) 

and 

1 7 P Pro P 

r p ~ r^ p ~ p1 0 ~ r 

At the downstream boundary, the static pressure is prescribed so that 

(2.16) 

P ~ Pexit 


(2.17) 

which when linearized becomes 



P = Pexit “ P 


(2.18) 


Although the problems considered in this report all have subsonic inflow and outflow, 
supersonic inflows and outflows could also be examined. For supersonic inflow, all three 
primitive variables would be specified at the inflow boundary. If the flow were supersonic 
at the exit, no boundary conditions would be specifed at the exit. 

2.2.4 Assembly of the Discrete Equations 

The two inlet and one exit boundary conditions, together with the 3 (/ — 1) con- 
servation equations, give 3 1 equations in 3 1 unknowns which completely specifies the 
problem. Now it is simply a matter of solving a set of linear equations to determine 
the perturbation variables. The most obvious way to set up the matrix equation to be 
solved is shown in Figure 2.2. The top two rows of the matrix equation are the inlet 
boundary conditions. Below them are the / — 1 sets of conservation equations, three for 
each cell. Finally, at the bottom of the matrix is the single exit boundary condition. 
This arrangement works quite well providing a small bandwidth matrix equation to be 
solved by Gaussian elimination. However, a slightly different approach will be used to 
assemble the matrix equations. 
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Figure 2.2: Conventional assembly of matrix equations for solution of the Newton iter- 
ation equations 
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Let the portion of the matrix having to do with the conservation equations is denoted 
by M and the corresponding right hand side by f . Similarly the portion of the matrix 
having to do with the boundary conditions be denoted by M bc aQ d the corresponding 
right hand side by f bc- Then one can formulate the solution of the linear equations as 
a minimization problem. That is to say, we wish to minimize the the sum of the squares 
of the residuals of the conservation equations plus the sum of the boundary condition 
residuals. This is equivalent to the statement that 

[M r M + M5 c M fl c]u = M T f + M5 c f B c (2.19) 


An Example 

To demonstrate the use of this technique, consider the one-dimensional problem 

§£ = /, *€[0,1], «(0) = U (2.20) 

This equation can be discretized using a central difference operator to give 

= A*/, + ^ (2.21) 

where 

A* = x i+ i - Xi 

/ 1+i = i (/< + /«■!) 

Suppose that the domain is composed of three cells (four nodes) so that in matrix form, 
Equation 2.21 can be expressed as 


-110 0 
0-110 
0 0-11 



«1 


- 


l»J 


Az/i + i 




«s 

. “4 . 



( 2 . 22 ) 


where A* = 1/3. Note that at this point, there are fewer equations than unknowns 
since the boundary conditions are not included. For the present example, the boundary 
condition is applied at the left boundary so that 


t n = U 


(2.23) 
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or in matrix form 


[ 1 0 0 0 


(2.24) 


The straightforward approach would be to place this equation above the the field equa- 


tions to give 


10 0 0 
■110 0 


(2.25) 


-110 0 t »2 Ax /,,1 

1+ * (2.25) 

0—110 uj Ax/ 2+ | 

0 0 *-l 1 J [ u 4 ] [ Ax/ 3+ i 

Using the least-squares approach, however, the matrix is formed by minimizing the sum 
of the square of the residuals of the difference equations and the boundary conditions. 
Hence, 


1 0 0 0 + 


-10 0 
1 -1 0 


-110 0 
0-110 


0 1 


-10 0 
1-10 


After simplification, this becomes 


2-100 

-12-10 

0-12-1 


0 0-11 


1 -1 


0 0 


Al /i+* 

Ax f 2+1 

A*/ 3+ i 


U - A x/ l+i 

A * A +* " 

Ax fi+$ ~ Ax ft+\ 


Note that the matrix on the left-hand side is positive definite. 
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Consider > for example, the case of / = 1 and U — 1. Then the exact solution to 
u* = 1 is u = x + 1. It is quite simple to verify that the above matrix equation gives 
the exact solution of u = [1, j, |,2]. 

Interpretation of the Squaring Technique 

The method described above resembles the least-squares finite element method [28,29]. 
In the finite element method, the solution to a differential equation is approximated 
by local shape functions. So for example, the approximate solution anywhere within 
a cell or element can be expressed in terms of the value of the function at the nodes 
of the element. The values of the discrete nodal values are then chosen such that the 
approximate solution satifies — in some sense — the partial differential equations. For 
the least-squares finite element method, the differential equation is “satisfied” if the 
residuals of the differential equation is minimized over the domain and the boundary 
condition residuals are minimized. In fact, the least-squares technique can be thought 
of as a variational principle for non-self-adjoint differential equations. 

Originally in this research, least-squares finite elements were used to discretize the 
PDEs which describe the steady and unsteady channel flow — with dissapointing results. 
We found, as did Astley, Walkington, and Eversman [30], that the method as described 
above gives poor results unless higher order shape functions are used. The present 
method more nearly resembles collocated least squares } where the residuals are minimized 
at a finite number of stations rather than over a continuous region. 

There are several advantages to the least-squares technique. Because the resulting 
equations are positive definite, it may be possible to make use of iterative matrix solvers 
which takes advantage of positive definiteness and sparseness, although in the present 
work only direct solvers were used. The best feature, however, is the ease with which 
boundary conditions are introduced. The assembly of the discrete field equations is 
independent of the boundary conditions. To introduce the boundary conditions, The 
squared boundary condition equations are simply added to the squared field equations. 
This is especially useful in the two-dimensional problem where reordering would be quite 
complicated. 
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2.2.5 Subsonic Flow in a Converging-Diverging Channel 

In this section, a model channel flow problem will be solved using the method out- 
lined above. The channel has an area distribution given by A(x) = 1 — jsin 2 (jrx) on 
the domain x = [0,1]. The flow is from left to right with the inlet toted pressure Pro 
and density pro specified at the inlet, and the static pressure p ex ^ specified at the exit. 
For this case the total pressure Pro and total density pro are 1.175 and 1.122. The exit 
pressure p ex j^ is 1.0. At this back pressure, the channel is not choked. A computational 
grid with 21 nodes was used (Ax = 0.05). Since the Euler equations are solved with a 
Newton iteration procedure, an initial solution must be assumed. The initial flow was 
assumed to be uniform with p = 1, U = 0.5, and P = 1. The flow was calculated using 
the Newton iteration procedure in four iterations. The converged solution is shown in 
Figure 2.3 along with the exact solution. Not surprisingly, the two solutions are nearly 
identical. The second-order accurate scheme allows accurate results to be computed 
with relatively few nodes. 

As is characteristic with the Newton iteration procedure, convergence is extremely 
fast when the solution estimate is not too far away from the actual solution. One draw- 
back, however, is that the method is not guaranteed to converge and in fact can diverge 
quite spectacularly if a poor initial solution is chosen. Figure 2.4 shows the strong 
convergence history exhibited in this case. Plotted is the logarithm of the maximum 
residual versus the iteration number. One can see that initially, the convergence is 
slightly faster than linear on a semi-logarithmic scale. After three iterations, no further 
convergence is seen which is due to round-off error associated with finite computer pre- 
cision. Figure 2.5 shows the initial pressure along with the calculated pressure after 
each of the first five iterations. This provides another illustration of the extremely fast 
convergence. Especially after the first iteration, the solution converges to its final value 
very quickly. 

2.2.0 Steady Shock Fitting 

In the previous section, we calculated the steady subsonic flow through a converging- 
diverging channel. If the flow is transonic, however, the presence of a sonic point and 
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Figure 2.3: Calculated steady flow through converging-diverging channel. Also shown 
for comparison is the exact solution. Total inlet pressure to exit pressure ratio is 1.175. 
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Figure 2.5: Computed pressure distribution history after each iteration. Note solution 
is converged after four iterations. 
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shock wave introduce additional complications to the problem. In this section, the 
numerical procedure for fitting the shock will be discussed, as well as the need for some 
smoothing to handle the sonic line. 

The Euler equations are nonlinear partial differential equations and as such they 
admit so-called weak solutions, i.e., solutions which contain jump discontinuities which 
satisfy the differential equations in an integral sense. The physical manifestation of 
this mathematical idealization is a shock wave. Consider once again the steady Euler 

equations. If these equations are satisfied everywhere, then 

fX.+t / Q d \ 

L. fe AP + A ai p ) rfx = 0 < 2 28 > 

where the integration runs from one side of the shock to the other. As e approaches 
zero, the limit of the integral gives the condition 

[[F + P]] = 0 (2.29) 

where the symbol [[• • •]] denotes the difference in the enclosed quantity across the dis- 
continuity. Equation 2.29 (the Rankine-Hugoniot equations) relates the flow on one side 
of the shock to that on the other side. 

These jump conditions can be linearized about a current estimate of the flow solution 
to obtain linearized shock-jump conditions for use with the Newton iteration procedure. 
These are given by 

[[(B, + B s )u|] = -[{F + P]] (2.30) 

and are applied at what will be the new estimate of the shock location. The term on the 
right hand side represents the jump condition at the new shock location calculated from 
the old flow estimate. This is not a particularly useful form of the jump equations since 
the new shock location is not known a priori. Instead, a new variable is introduced 
which is the distance the shock moves from the current shock location to the new 
shock location. With the aid of this additional perturbation variable, the shock jump 
conditions can be expressed in terms of the flow at the current shock position without a 
priori knowledge of the new shock position by extrapolating the shock jump conditions 
from the current shock position to the new shock position, i.e., 

|[(B, + B s ) u]] + ~ ([[(B, + B,) u]|) *. = - [|F + P)] - ^ ([(F + P])) *. (2.31) 
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where now this condition is applied at the current shock position rather than the new 
shock position. Note that the second term on the left-hand side of Equation 2.31 is 
second order and can be neglected. Hence retaining only first-order effects gives 

[[(Bj + B s )n+i.jj(F + P)]] =-||P + P]] (2.32) 

The numerical procedure is as follows. An initial shock position is assumed. The 
grid is then constructed as before except that at the shock there are two grid points, 
one grid point corresponding to the flow on the upstream side of the shock, and one 
corresponding to the flow on the downstream side. These two nodes comprise the shock 
cell. The conservation equations are applied to the conservation cells as before, but at 
the shock cell the linearized shock-jump conditions are applied with one-sided, three- 
point, second-order accurate differences used to evaluate the gradients in the base flow. 

Assuming the flow is subsonic at the inlet and exit, two upstream and one down- 
stream boundary conditions are imposed as in the subsonic example. But this does 
not completely specify the problem. Taking count of the equations and unknowns, we 
have 7 grid points each with three perturbation variables plus the extra shock position 
perturbation variable for a total of 37 + 1 unknowns. There are 3(7 - 2) conservation 
equations arising from the 7-2 conservation cells, three jump conditions at the shock 
cell, and two upstream and one downstream boundary condition for a total of 37 equa- 
tions. Unfortunately, we are one equation short. This is not due to the addition of 
the shock perturbation variable but rather to the presence of a sonic point in the flow. 
Giles [31] solved this problem by adding a special equation at the sonic point which 
describes the behavior of the J~ characteristic. A more common approach, which was 
also used in this research, is to add a small amount of smoothing to the equations to 
stabilize the sonic point. The smoothing is added to the matrix equation after the squar- 
ing procedure described earlier. A small amount of a Laplacian like operator is added 
to the "stiffness” matrix. The operator is positive semi-definite which preserves the 
symmetry of the equations and makes the system positive definite. The final equation 
matrix equations is 

[M T M + Ml c M B c + «L] u = M T f + Ml c f BC + (2.33) 
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where 


f L = -LU 

and e is a small parameter used to control the level of smoothing introduced into the 
scheme. 

The Newton iteration procedure is then as follows. An initial transonic flow is as- 
sumed, including the initial position of the shock. The linearized conservation equations, 
shock jump equations, and boundary conditions are set up and “squared” as described 
above. A small amount of Laplacian smoothing is then added to the system. The 
system of equations is solved using a sparse Gaussian elimination solver to obtain the 
perturbation primitive variables and the perturbation in the shock position. The prim- 
itive variables are updated, then the grid is readjusted so that the double-noded shock 
cell is coincident with the new shock position. This requires that the updated primitive 
variables be interpolated from the old grid to the new grid. The entire procedure is 
repeated until convergence is achieved. The procedure exhibits rapid convergence when 
it works, but convergence is not guaranteed. In practice the method usually converges 
in less than ten iterations if a reasonable initial guess is made. 

2.2.7 Steady Transonic Channel Flow 

To demonstrate the shock fitting algorithm, we will consider the flow through the 
same convergent divergent channel previously considered in Section 2.2.5. This time, 
however, the total pressure and density are somewhat higher so that the flow is choked. 
The total inlet pressure Fro and total density pro are 1.25 and 1.173, respectively. The 
exit static pressure p ex jj is 1.0. The initial shock position is assumed to be X, = 0.8. 
Recall that a poor initial guess can cause the procedure to diverge. Therefore, one 
should use whatever experience and knowledge one has about the nature of the actual 
solution to pick as reasonable a first guess as possible. The initial starting solution 
for the Newton iteration procedure is shown in Figure 2.6. Starting with this initial 
solution, the Newton iteration procedure was run for enough iterations to insure a 
converged solution. The converged solution is shown in Figure 2.7 along with the exact 
solution. Again excellent agreement is seen between the computed solution and the 
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Figure 2.7: Converged solution for choked flow in one-dimensional channel. Total inlet 
to exit pressure ration is 1.25. 
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Figure 2.8: Convergence history of choked flow calculation. Initially the solution “wan- 
ders” but then converges very quickly. 

exact solution. Note in particular the precise definition of the shock arising from the 
use of shock fitting. The shock position is calculated to be X t = 0.6805 while the exact 
shock position is X, = 0.6788. This accurate an estimate of the shock position could not 
be obtained using shock capturing with the same grid resolution since the uncertainty 
in a captured shock’s position is on the order of a few conservation cells. 

As with the subsonic steady flow calculation, convergence is extremely fast. The 
maximum residuals after each iteration are shown in Figure 2.8. Note that the con- 
vergence history can be divided into roughly three phases. In the first phase, the 
convergence is slow and erratic, with the residuals even increasing on occasion. In the 
second phase, the residuals decrease rapidly and nearly linearly on a logarithmic scale. 
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Beyond this point, round off errors preclude any improvement in the solution and the 
computation is for all intents converged. The three phases of convergence can also be 
seen in Figure 2.9. Plotted is the computed shock position after each iteration as well 
as the perturbation in shock position. For the first four iterations, the shock position 
“bounces” around quite a bit. This corresponds to the initial phase of convergence. If 
the computation is to diverge, it will usually happen in the first few iterations. Once 
the solution gets in the neighborhood of the actual solution, the shock position quickly 
converges to its final value. After about ten iterations, the shock has settled in its final 
position, save for extremely small excursions due to round-off error. 

Because the numerical algorithm is second order accurate, the computed solution 
should differ from the exact solution by 0(Ax J ). As a test of this hypothesis, the grid 
resolution was varied from three conservation cells up to 96 conservation cells. For 
each case, the steady shock position was calculated. Figure 2.10 shows the error in 
the computed shock position as a function of the average cell size. Note that on a 
logarithmic scale, the error curve has a slope of about two, demonstrating the scheme is 
second order accurate. In particular, for a A* of about 0.04, the error in the computed 
shock position is less than 0.1 percent of the channel length. 

2.3 Solution of the Linearized Unsteady Euler Equations 

In the previous section, steady subsonic and steady transonic flows through a one- 
dimensional channel were calculated. In this section, unsteady perturbation flows about 
these mean flows will be calculated. These solutions will be compared against a time- 
accurate time-marching scheme to verify the results. Although the flows calculated 
here are fairly simple one-dimensional flows, they will demonstrate many of the salient 
features of the present analysis. 
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Figure 2.9: Computed shock position history. The upper plot shows the shock position 
after each iteration. The lower plot shows change in position from one iteration to the 





2.3.1 Discretization of the Linearized Euler Equations 

The linearized Euler equations (Equation 2.7) are discretized using a cell centered 
finite difference scheme given by 

. Ui+I+Ui , (AB : u) t+ i - (ABju)^ 

B,) j+ j + ^ + 

A, +i (BjU) '^ (BlU) ' = 0 (2.34) 

These difference equation are second order accurate in Ax. Note also that they are 
complex due to the harmonic assumption. It can be shown, at least in the constant 
area case, that these equations are nondissipative. That is to say that the pressure 
and entropy waves do not grow or decay. This is due to the central differencing. The 
equations are, however, dispersive. The waves may not “travel" at quite the right speed. 
This is especially true of the shorter wavelength waves. 

2.3.2 Unsteady Channel Flow Boundary Conditions 

In the two-dimensional cascade flow problem, the unsteady boundary conditions re- 
flect will real physical problems encountered in turbomachinery such as inlet distortion, 
potential disturbances from neighboring blade rows, incoming vorticity and entropy 
waves, or motion of the cascade blades themselves. In these one-dimensional model 
problems, the boundary conditions are somewhat contrived, but must be mathemat- 
ically well posed. If the flow is subsonic at the inlet and exit, then two boundary 
conditions will be required at the inlet, and one will be required at the exit. If the flow 
is supersonic, three are required at the inlet and none are required at the exit. 

In the examples considered in this chapter, the flow at he inlet and exit are subsonic, 

and hence, the unsteady total pressure and total density will be specified at the inlet 

and the unsteady static pressure will be specified at the exit. At the inlet the boundary 

conditions can be expressed as 

7Pr _ Pro + _ 7P t 7 - 1 -2 

PT ^ Pto + A pT*i ut P 2 

Pr Pr 0 -f- A/Ve ,wt P 

Pt ( PTO + A PT« ,h ’*) 1 P 1 


(2.35) 

(2.36) 
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which when linearized become 


and 


p’- r +1# ' pJoV/Vo 


A Pp A pt\ 

Pro pro / 

(2.37) 

7Ap r \ 
Pro / 

(2.38) 


jP-fiTTP 

Note in particular that if the inlet low is isentropic (to first order), the expression on 
the right hand side of Equation 2.38 vanishes. Similarly, the downstream static pressure 
boundary condition is 

P = Pexit = ^exit + (2.39) 

which when linearized becomes 

p = AP exit (2.40) 


2.3.3 Assembly of the Discrete Equations 

The numerical boundary conditions together with the finite difference equations 
completely specify the linear equations to be solved. As with the Newton iteration 
equations, the equations are “squared” to produce a positive definite matrix equation 
given by 

[M H M + Mf c M BC ] u = Mg c f BC (2.41) 

where the transpose operator has been replaced by the Hermitian operator because the 
original matrices are complex. The resulting positive definite matrix is Hermitian. 


2.3.4 Unsteady Subsonic Results 

As an example, we will compute the small disturbance unsteady flow in a channel. 
The steady flow about which the Euler equations are linearized is the subsonic flow 
considered in Section 2.2.5. Recall that the inlet total pressure Pro is 1.175, and the 
inlet total density pro is 1.122. Downstream the static pressure P ex j t is 1.0. Because 
the linearized unsteady Euler equations are linear and homogeneous, the calculated flow 
scales with the unsteady boundary conditions. For convenience then, The total pressure 
at the inlet varies as APr = 1-0. The excitation frequency w is 1.257. Since the steady 
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Figure 2.11: Unsteady pressure distribution due to upstream total pressure and density 
variation. Excitation frequency w is 1.047. 

velocity in the channel and the channel length are order unity, one can think of this 
frequency as a reduced frequency. The total density is varied such that the inlet flow is 
isentropic (A Pr/Pr = 1&Pt/pto)- Downstream the static pressure variation AP ex j t is 
set to zero (this is a reflecting boundary condition). The calculated unsteady pressure 
distribution is shown in Figure 2.11. Because the perturbation variables are complex, 
both the real and imaginary parts of the pressure must be presented to completely 
describe the pressure distribution. Alternatively, one could present the magnitude and 
phase of the pressure at each point. The former is preferred, however, since physically 
the real and ima gin ary parts each represent a snapshot of the flow at a particular instant 
in time. 

To validate the linearized Euler method, these results were compared to those ob- 
tained from the time-accurate, time-marching code of Giles [31]. This code is unique in 
that shocks are fit rather than captured, and no artificial viscosity is required although 
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special treatment is required at sonic points. The procedure was to first time march 
the code with steady state boundary conditions until the solution had nearly converged. 
Then the code was restarted, but this time the upstream total pressure was varied si- 
nusoidally in time with the total density chosen such that the inlet flow was isentropic. 
Downstream, the static pressure was held constant. Because the time-marching code is 
nonlinear, the total pressure variation was chosen to be very small so that the results 
can be compared to linearized theory. (In the next chapter, the effect nonlinearities due 
to larger levels of unsteadiness will be examined.) After any initial unsteady transients 
have decayed, the unsteady flow is periodic and very nearly sinusoidal. The time history 
was then Fourier transformed to convert the time domain results into the frequency do- 
main so that they can be compared directly to the linearized Euler results. These time 
marching results are compared to the linearized Euler results in Figure 2.11. Note the 
excellent agreement between the two methods. 

For the next case, the frequency w is increased to 2.094 The calculated pressure 
distribution is shown in Figure 2.12. Again, good agreement is seen between the 
linearized Euler theory and the time-marching code. 

These results demonstrate the validity of the small disturbance approach. If one is 
interested in small disturbance flows which are nearly harmonic in time, the linearized 
Euler method offers an alternative to the time marching Euler approach. Furthermore, 
the linearized Euler method requires less computational time than the time marching 
technique. In one dimension, the CPU time is reduced by three orders of magnitude, 
while in two dimensions, the CPU time reduction is more moderate at one to two orders 
of magnitude. On the other hand, the linearized Euler method is not suitable for flows 
with large disturbances, or flows which are not periodic. 

2.3.5 Unsteady Shock Fitting 

If the mean flow is transonic, then the calculation of the unsteady flow has the added 
complexity of a moving shock. The complexity arises because the shock is moving, but 
the grid is stationary. Therefore, just as in the case of the Newton iteration procedure 
for shock fitting, we will need some form of extrapolation to accurately model the shock 
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Figure 2.12: Unsteady pressure distribution due to upstream total pressure and density 
variation. Excitation frequency oj is 2.094. 
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motion. 


The starting point of the analysis are the full nonlinear unsteady shock jump con- 
ditions. To derive them, a change in coordinate systems is made. The new coordinate 
system is fixed to the moving shock. Hence 

x = x - x,(t) (2-42) 

t = t (2.43) 

using the chain rule of differentiation yields 


8*1* 

II 

£1* 

(2.44) 

a a a 

dt dt V * dx 

(2.45) 


where v, is the instantaneous shock speed in the fixed coordinate system. Using the 
above expressions for the differential operators to derive the unsteady Euler equations 
in the shock coordinate system gives 

(a-'») AU+ w" +A 5s p -° (246) 

Integration of this equation from just upstream to just downstream of the shock gives 
the desired shock jump conditions. Because the time derivative is bounded in the shock 
frame of reference and the interval of the integration goes to zero, the time derivative 
does not appear in the shock-jump condition. The spatial derivatives, however, act 
like impulses at the shock, and hence, integration of these impulses yields the jump 
condition 

[[F + P-v,U]] = 0 (2.47) 

These jump conditions give the nonlinear relationship between the flow on either side 
of a moving shock. These are the unsteady Rankine-Hugoniot equations. 

Next, let us assume that like the primitive variables, the shock motion is also har- 
monic in time. Hence the shock position is composed of mean part plus a harmonic 
variation. 

x.{t) = X, + x,e’ ut (2.48) 
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Therefore, the instantaneous shock velocity is x,e jut . Linearizing the shock jump con- 
dition gives 

[[(B 2 + B s ) u - jojx,V ]] = 0 (2.49) 

It is emphasized that the jump condition is applied at the instantaneous position of the 
shock. As with the Newton iteration shock jump equations, we can extrapolate from 
the nominal shock position to the instantaneous shock position and apply the shock 
jump condition at the nominal shock. Keeping only first order terms gives 

[[(B, + B s ) u + x.^ (F + P) - jwz.u]] = 0 (2.50) 

2.3.6 Assembly of Unsteady Transonic Equations 

The matrix equations are assembled in essentially the same way the steady transonic 
equations were assembled. Note that because of the transition form subsonic to super- 
sonic flow at the sonic line, there will not be quite enough equations for unknowns. 
Hence, a small amount of smoothing is introduced to stabilize the sonic point. The 
resulting matrix equations are given by 

[m h M + + *1*] u = (2-51) 

2.3.7 Unsteady Transonic Results 

As an example, the unsteady small disturbance flow about a steady transonic flow 
through a converging diverging channel is considered. The base flow is the one consid- 
ered in Section 2.2.7. The inlet total pressure Pro is 1.25 while the inlet total density pro 
is 1.173. Downstream the static pressure P exit is 1.0. For the perturbation flow, the un- 
steady variation in the upstream total pressure A Pr is 1.0, with a frequency w of 1.257. 
The total density Apr is chosen such that the inlet flow is isentropic. Downstream, the 
variation in static pressure AP ex j^ is set to zero. The computed unsteady pressure distri- 
bution is shown in Figure 2.13. Also plotted is the time-marching solution of Giles [31]. 
Note the generally good agreement between the two theories except at the throat the 
time-marching code has a localized problem. Nevertheless, away from the throat, the 
agreement is excellent. The shock motion is also in good agreement. The complex shock 
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Figure 2.13: Unsteady pressure distribution for choked flow in one-dimensional channel. 
Upstream total density and pressure is varied with a frequency u> = 1.257 
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Figure 2.14: Unsteady pressure distribution for a downstream pressure perturbation 
with frequency oj of 3.1416. Note that ahead of the shock, the pressure is zero. 

motion as computed by the linearized Euler theory is x $ = —.1724 — 0.8522j while the 
equivalent shock motion predicted by the time-marching code is x 9 = —.1707 — 0.8697j. 

As a final example, consider the case of a downstream pressure perturbation. Up- 
stream the total pressure and density perturbations are set to zero. The excitation 
frequency u) is 3.1416. The computed pressure distribution is shown in Figure 2.14. 
Note again the excellent agreement between the linearized Euler theory and the time- 
marching Euler theory. Also, ahead of the shock, the pressure is zero for both theories 
as it should be, since small disturbances cannot travel upstream in a supersonic region. 
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2.4 Summary 


In this chapter, the numerical method for calculating small disturbance flows in a 
one-dimensional channel was presented. Although only one-dimensional model problems 
were analyzed, the theory is essentially the same for two-dimensional problems which 
we will consider in later chapters. 

Also presented was a Newton iteration procedure for solving for the steady channel 
flow. The method is extremely efficient with subsonic problems typically converging in 
about five iterations, and transonic problems converging in about ten iterations. Fur- 
thermore, since the Newton iteration procedure is so similar to the unsteady linearized 
Euler method, there is a significant savings in code development time. The steady code 
is second-order accurate and computed steady flows agree well with exact solutions. 

Unsteady pressure distributions were calculated for several model unsteady flow 
problems. These solutions were compared to a time-accurate time-marching scheme 
with excellent agreement. 

Shock fitting was used to model the position of both steady and unsteady shocks. 
This allows for a sharp definition of the shock with relatively few conservation cells. 
Furthermore, unlike captured shocks, the computed shock position for the steady flow 
and the computed shock motion for the unsteady flow are second order accurate. In 
the steady case, the computed shock position is nearly indistinguishable from the exact 
shock position. In the unsteady case, the unsteady shock motion agrees well with the 
unsteady motion computed using a time-accurate time-marching scheme. 

In the next chapter, the differences between the linearized Euler theory and the lin- 
earized full potential theory will be examined, as well as the limitations the linearization 
places on the linearized Euler theory. 
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Chapter 3 


Advantages and Limitations of 
the Linearized Euler Analysis 


In the previous chapter, the linearized Euler anlaysis used to predict small dis- 
turbance, quasi-one-dimensional unsteady flows was presented. In later chapters, the 
method will be extended to two-dimensional channel and cascade flows. Before pro- 
ceeding, however, two important issues will be addressed in this chapter. The first issue 
is whether the linearized Euler method is useful for evaluating unsteady flow problems 
with finite levels of unsteadiness in the flow. The second is whether the method does 
a significantly better job at predicting unsteady flows than the currently available and 
less computationally expensive full potential methods. 

The fundamental assumption on which the linearized analyses are based is that the 
flow can be divided into two parts: a mean flow component and a much smaller unsteady 
component. The linearized analysis asymptotically approaches the full time-dependent 
nonlinear analysis as the unsteadiness goes to zero. But if the unsteadiness is finite, 
nonlinear effects may also be important. In Section 3.1, we will examine the behavior 
of the full nonlinear Euler equations for various levels of unsteadiness. An asymptotic 
analysis of the Euler equations shows that the linear assumption is valid unless the 
square of the unsteadiness is significant. Numerical examples will demonstrate that the 
linearized Euler analysis does a good job of predicting the unsteady behavior even for 


49 



moderate levels of unsteadiness. 

The second issue is whether the linearized Euler analysis provides a significant im- 
provement in modelling of unsteady flows over the current generation of linearized un- 
steady full potential methods. Because of the larger computational effort required to 
perform the linearized Euler analysis, its use would be unwarranted for cases in which 
the potential method gives accurate predictions. For flows which are isentropic and 
irrotational, the two methods in principle will provide identical predictions. For non- 
isentropic flows, such as flows with shocks, the two methods will not agree. The question 
is how strong the shocks must be before the two methods begin to differ significantly, 
and what are the mechanisms which produce these differences. In Section 3.2, we will 
show that the isentropic assumption inherent in the potential methods can produce er- 
rors which are significant even for flows with weak shocks. Furthermore, the potential 
methods are shown to become unreliable for shocked flows at low reduced frequencies. 


3.1 Nonlinearities in Unsteady Flow 

Much of the work in unsteady flows in cascades has been based on the assumption 
that the fluid flow in a cascade can be linearized about some nominal steady flow. For 
relatively small perturbations in the flow, this is a reasonable assumption and certainly 
goes a long way toward simplifying the nonlinear unsteady compressible flow equations. 
But the perturbation analysis may not be valid for flows with finite perturbations or for 
the levels of unsteadiness typically found in turbomachines. In this section, the effect of 
nonlinearities will be examined using two approaches. In Section 3.1.1, a Fourier series 
analysis is used to demonstrate the appearence of higher harmonic components due 
to nonlinearities. In Section 3.1.2, a series of numerical experiments will be performed 
using a time-accurate time-marching Euler code to quantitatively demonstrate the effect 
of nonlinearities on unsteady flow solutions for finite flow perturbations. 
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3.1.1 Fourier Analysis of Unsteady Flows 

To begin our discussion, consider again the continuity equation for flow in a quasi- 
one-dimensional channel. 

dp d 

A m + Tt puA = 0 < 31 > 

The continuity equation, like the momentum and energy equations, is nonlinear in the 
primitive variables and homogeneous. Since the Euler equations are homogeneous, any 
unsteadiness will arise from the applied boundary conditions such as those due to up- 
stream flow perturbations or downstream pressure perturbations. Suppose that these 
external excitations are sinusoidal in time with frequency w. Because the Euler equa- 
tions are nonlinear, the exact time dependent solution in general will not be sinusoidal. 
However, after any initial transients which arise due to the initial conditions decay 
away, the flow will be periodic in time with period T = 2ir /u>. Therefore, the primi- 
tive variables can be expanded in a complex Fourier series with fundamental frequency 
w [32]. 

/>(*,«)= £ />«(*)«*"* (3.2) 

n=-oo 

«(!,<)= £ «„(*)«*•"■• (3.3) 

n=— oo 

P{x,t) = £ Pn{x)e’ unt (3.4) 

n=-oo 

where p n , u n , and p„ are complex functions of * alone. 

To obtain the Fourier coefficients, a Fourier transform is used to convert the time 
history of the periodic flow into the frequency domain. For complex Fourier series, the 
transform is given by 

Po(*) = fJ Q p{ x ^) dt (3 5) 

Pn(x) = f p{x,t)e~’ unt dt (3.6) 

Alternatively, one could represent the periodic time history of the flow in a Fourier 
sine and cosine series. Hence, for example, the density is given by 

OO OO 

p[x, t) = p 0 {x) + Yi Pen[x ) cos(nu/t) + Y P*n{x) sin(nwt) (3.7) 

n=l n=l 
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where the coefficients p 0 , Pen , and Pm are given by 


Po{x) - 


(3.8) 

Pcn{x) = | j 

rT 

1 p(x, t) cos(wnt) dt 
'o 

(3.9) 

Psn(x) = ji ' 

rT 

1 p(x,t)sin(unt) dt 
lo 

(3.10) 


The Fourier sine and cosine series is really equivalent to the complex Fourier series. 
The relationship between the complex Fourier series coefficients and the sine and cosine 
series coefficients is 


Pen{x) = (p n {x) + P-n{x)) 

(3.11) 

Pm(x) = 7 ( Pn(x ) - P-n(x)) 

3 

(3.12) 

Furthermore, if the function p(x,t) is a real function, then 


Pen{x) — R^Pnt*)) 

(3.13) 

P»n{x) = -Im(p n (*)) 

(3.14) 

P-n = Pn 

(3.15) 


where the * symbol indicates the conjugate operator. 

So long as the flow is periodic, the Fourier series representation is a perfectly general 
way to describe the flow. More importantly, however, is the fact that the Fourier series 
representation allows us to systematically examine the effects of nonlinearities in the 
flow. The general procedure is to substitute the complex Fourier series representation 
of the primitive variables into the Euler equations. The resulting partial differential 
equations are then in terms of the Fourier coefficients. Hence, the differential equations 
describe the way the coefficients are related to one another. A careful asymptotic 
analysis of the behavior of the Fourier coefficients provides insight into the behavior of 
the nonlinear flow. 

To begin, Equations 3.2-3.4 are substituted into the Euler equations. Wherever 
products of primitive variables occur, the multiplication of the Fourier series is carried 
out term by term. For example, in the continuity equation the multiplication of the 
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product pu is carried out term by term and the resulting series is rearranged in like 
powers of exp(jwt). For the continuity equation to hold, the resulting series should 
vanish term by term. Hence we have that 

Q 

—A {pqxjlo + Piti-i + P 2 U -2 + • • • + P- i«i + P-2 U 2 + •••) = 0 (3.16) 

ox 

Q 

juApi + —A (po w l + Pi^o + P2 u -i + • • ■ + P-1^2 + P- 2 U 3 + •••) = 0 (3*17) 

ox 

Q 

— juAp-i + — A (po u -i + P- l^o + P- 2 «i + • • • + Pi«-j + P 2 «-s H ) = 0 (3.18) 

ox 

d 

2 jwAp 2 + —A (piuj + P2U0 + P3W-1 + • • • + P0U2 + P-1^3 + * • •) = 0 ( 3 . 19 ) 

ox 

Q 

— 2jujAp-i + —A (p_iu_i + p_ 2«o + P-s^i + • • • + P 0 V -2 + P 1 U -3 + • • •) = 0 (3.20) 

ox 

and so forth. 

To solve completely for the periodic flow in the channel, one would have to solve 
simultaneously for all the Fourier coefficients. This of course is impossible since in 
principle one would have to solve an infinite number of nonlinear simultaneous differ- 
ential equations. Nevertheless, some insight into the effects of nonlinearities on higher 
harmonic content can be found by looking at the asymptotic behavior of the Fourier 
components as the perturbation level gets small. Let us assume for the moment that 
Po » P±i P±n and uo 3> u±i » u± n . Then to leading order, the Fourier coefficients 
are given by 

J^pouoA = 0 (3.21) 

d 

juApi + — A(poUi + piuo) = 0 (3.22) 

Q 

2juAp2 + —A (piU! + p 2 «o + P0«2) = 0 (3.23) 

Equation 3.21 is simply the steady state continuity equation. That is to say that the 
leading order behavior of the mean flow continuity equation is equivalent to the steady 
flow continuity equation. Equation 3.22 is the linearized continuity equation on which 
the unsteady analysis developed in this report is based. Finally, Equation 3.23 is the 
leading order behavior of the second harmonic component in the flow. Note that the 
second harmonic component is excited by the perturbation product of the fundamental 
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components p%ui. This means that for small levels of unsteady flow, the unsteady 
motion of the fluid will be nearly harmonic in time as this product is extremely small. 
However, as the level of unsteadiness gets larger, the product will become larger until at 
some point the second harmonic content will become significant. The component of the 
second fundamental frequency present will be approximately proportional to the square 
of the excitation level. 

We have just seen that to leading order, the linearized Euler equations give the cor* 
rect results. Next consider the higher order corrections to Equations 3.21-3.23. Keeping 
only the next highest order terms in the mean flow equation gives 

— ^Ap 0 UQ + ~(p e i« c i + = 0 (3.24) 

In other words, the mean flow is altered by the presence of unsteadiness. This effect, 
although second order, can be important if the unsteadiness is large. These modifying 
terms in the continuity equation (and also in the momentum and energy equations) are 
analogous to the Reynolds stress terms which appear in the analysis of turbulent flows. 
Because the equations which describe the unsteady flow are nonlinear, the mean flow is 
not the same as the steady flow. For boundary layers, this means that the mean flow for 
a laminar (steady) flow is not the same as the mean flow for a turbulent (unsteady) flow. 
For aeroelastic problems, it means that the steady pressure distribution on an airfoil is 
not the same as the mean pressure distribution on a vibrating airfoil. Fortunately, this 
effect is second order and therefore vanishes as the unsteadiness becomes small. 

Similarly, the higher order corrections to the first order (linearized) Euler equations 
are third order, or two orders higher than the leading order behavior. Therefore the 
linearized equations are valid so long as the square of the unsteadiness is not significant. 
Said another way, the difference between the results of the linearized Euler analysis and 
an exact solution will be of the order of the square of the unsteadiness. This sort of 
asymptotic analysis shows qualitatively that the linearized Euler analysis is valid for 
small to moderate levels of unsteadiness. In the next section, a time-marching Euler 
code will be used to determine the unsteady flow in a channel for increasing levels of 
unsteadiness. The results will be compared to the qualitative analysis of this section. 


54 



3.1.2 Numerical Examples 


In Section 3.1.2, the qualitative behavior of nonlinearities in unsteady flow was 
examined by looking at the asymptotic behavior of the various Fourier components 
of the flow. In this section, we will numerically examine the effects of nonlinearities 
on unsteady flows by performing a numerical experiment. In this experiment, a time- 
accurate Euler code is used to calculate the unsteady flow in a quasi-one-dimensional 
channel. The source of the unsteadiness in the flow is a harmonic variation in the 
upstream total pressure and density. The calculated time history of the flow can then 
be Fourier transformed to obtain the Fourier components of the flow. 

For this numerical experiment, the one-dimensional, time-accurate, conservative, 
implicit Euler code of Giles [31] was used as the test vehicle. The general procedure 
is to first run the code with steady boundary conditions at the inlet and exit of the 
channel. The code is time marched until the flow reaches a steady state solution. Then, 
using this solution as a starting point, the code is restarted, but this time the boundary 
conditions are sinusoidal in time. The code is time marched sufficiently long enough so 
that initial transients will have decayed. The resulting unsteady flow will be periodic 
in time although not necessarily sinusoidal due to the nonlinear nature of the Euler 
equations. The time history of one period of the flow is then Fourier transformed to 
obtain the zeroth, first, and second Fourier components of the flow. 

For the first case studied, the nominal mean flow was subsonic. The channel is 
defined on the domain x G [0,1] and has the area distribution A(x) = 1 — ^sin 2 (xx). 
Shown in Figure 2.3 are the nominal steady pressure, density, and velocity through the 
channel. For this case, the upstream total pressure Ppo was set to 1.175, the upstream 
total density pro to 1.122, and the downstream static pressure F e *it to unity. Giles’ 
code was time marched until a steady state solution was obtained. 

Next, to investigate the unsteady flow in the channel, the code was restarted and 
the upstream total pressure upstream was varied as Pr = Pro + A Pp cos(u >t) with a 
frequency u / of 1.0472. The upstream total density was varied such that the inlet flow 
was isentropic. For the cases reported here, the total pressure variation A Pp was 0.005, 
0.01, 0.02, 0.05, and 0.10. Downstream, the pressure was held fixed at unity. Giles’ 
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Figure 3.1: Pressure at the throat of a converging-diverging channel for an upstream 
variation in total pressure and density. The excitation frequency w is 1.0472. Note the 
higher harmonics present in the flow. 

code was time marched long enough for the flow to reach a periodic state, that is, until 
all transient behavior dependent on the initial conditions had decayed away. Then, one 
period of the response was Fourier transformed to determine the harmonic content of 
the flow. 

Figure 3.1 shows a typical calculated time history of the flow in the channel. For this 
case, APr = 0.05, and w = 1.0472. Plotted is the calculated pressure as a function of 
time at the throat of the channel. Clearly, the flow is not sinusoidal in time, but rather 
is periodic (or nearly so after the transients have passed) and contains a significant level 
of higher harmonic content. Figure 3.2 shows the computed harmonic content of the 
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Figure 3.2: Harmonic content of unsteady subsonic flow in a converging-diverging chan- 
nel for an upstream variation in total pressure and density. The excitation frequency w 
is 1.0472. Note the significant levels of the second harmonic. 
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Figure 3.3: Mean flow through a converging-diverging channel for various levels of 
upstream excitation. The excitation frequency u> is 1.0472 

zeroth, first, and second harmonic frequencies. The large level of the higher harmonic 
(second) content indicates that for this level of unsteadiness in the flow, nonlinear effects 
are important, at least if one wishes to determine the precise time history of the flow. 

Figure 3.3 shows the mean flow through the channel for various levels of upstream 
total pressure variation. Note that for total pressure variations of 0.005, 0.010, and 
0.020, the mean pressures are nearly indistinguishable from one another. For a total 
pressure variation of 0.050, the mean pressure deviates slightly from the steady solution 
at the throat. For a toted pressure variation of 0.100, the difference is even more pro- 
nounced. Recall that the mean flow is modified by a Reynolds stress like term in regions 
of large flow unsteadiness. As will be seen shortly, the region of greatest unsteadiness 
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Figure 3.4: First harmonic component in unsteady pressure for increasing levels of 
unsteadiness in the flow. 

is in the neighborhood of the throat. Hence, this is where nonlinear effects will first be 
noticed. 

Shown in Figure 3.4 are the first harmonic components calculated from Giles’ results. 
So that they all appear on the same scale, the results are scaled by the excitation level 
APr- The linear theory and the time-marching solution for a total pressure variation 
APr of 0.005, 0.01, and 0.02 give virtually identical results. Even the A Pp = 0.05 
case, which has significant levels of higher (second) harmonic components, has good 
agreement with the linearized theory. For higher levels of unsteadiness, however, the 
results begin to differ significantly. 

Next we wish to test the hypothesis that the second harmonic components we pro- 
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portional to the square of the excitation level for low levels of unsteadiness. Figure 3.5 
shows the second harmonic components scaled by A Pj>. Note that for APr = 0.005, 
0.01, and 0.02, the components agree well with one another. For higher levels of un- 
steadiness, the agreement is not as good, but are still of the same order. These results 
demonstrate clearly that for small levels of unsteadiness, the second harmonic terms 
appear with the square of the level of unsteadiness. 

Note that in all three cases, the zeroth, first and second harmonic components begin 
to deviate from the small their limiting values at a total pressure variation APr of about 
0.05. At this level of inlet disturbance, the magnitude of the first harmonic component 
of the pressure is a maximum near the throat and has a magnitude of about 0.075 or 

m 

about ten percent of the local pressure. The second harmonic component is around four 
percent of the steady pressure at the throat. Although one case is hardly representative 
of all unsteady flows, this example does give a feeling for the level of unsteadiness at 
which the linearized Euler analysis is no longer valid. 

As a second example, consider the transonic channel flow shown in Figure 2.7. The 
inlet total pressure and density are 1.25 and 1.173 respectively. The exit static pressure 
is 1.0. This back pressure is low enough to choke the flow. The flow enters the channel 
subsonically, passes through the sonic line at the throat, becomes supersonic, then 
shocks and exits the channel subsonically. Unsteadiness is introduced in this example 
by varying the back pressure sinusoidally in time such that p e »t = ^exit + Ap ex it cos(wt) 
where the frequency w is 1.0472. This unsteady back pressure produces unsteady flow 
aft of the shock as well as motion of the shock itself. The flow upstream of the shock 
will be quiescent unless the shock motion is so large that the shock moves upstream 
of the throat. The linear theory predicts that the unsteady pressure and shock motion 
will be sinusoidal. 

As with the subsonic case, Giles’ code was time marched until a steady-state solution 
was obtained. Then the code was restarted and run for various levels of unsteady back 
pressure. For the cases studied here, the unsteady variations of the back pressure Ap tx it 
are 0.010, 0.020, 0.040, and 0.080. For all cases except the last, the shock oscillated back 
and forth in the channel without passing through the throat. For the last case, however, 
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Figure 3.5: Second harmonic component in unsteady pressure for increasing levels of 
unsteadiness in the flow. 


61 




Table 3.1: Fourier components of nonlinear shock motion for various levels of unsteady 
back pressure. 


Apexit 

x, 0 

Re(*,i) 

Apexit 

Im(x,i) 

Apexit 

Re(* #a ) 

ApL* 

Im(», a ) 

A P»it 

0.01 

0.6791 

-0.0999 

1.1846 

0.9950 

-1.3732 

0.02 

0.6782 

-0.0934 

1.1918 

0.7949 

-1.8210 

0.04 

0.6737 

-0.0554 

1.2214 

0.7171 

-2.0072 

Linear 

0.6805 

-0.0996 

1.1777 

— 

— 


the shock motion was great enough so that the shock passed completely through the 
throat and then disappeared slightly ahead of the throat. At a later time in the cycle 
the shock reformed in the aft section of the channel and the cycle was repeated. 

The unsteady pressure time histories were Fourier transformed to obtain the zeroth-, 
first-, and second-harmonic components of the flow. These results for the zeroth- and 
first-harmonic are shown in Figures 3.6-3.7. The pressure in this region is anharmonic 
due the passing of the shock and the corresponding sudden pressure rise or fall. There- 
fore, the Fourier transform of this time history will resemble an impulse in the vicinity 
of the shock. The harmonic components of the time history are not presented in the 
region of the channel in which the shock oscillated. Figure 3.6 shows the mean flow 
through the channel along with the computed steady flow. Note that the mean flow is 
indistinguishable from the steady flow in all cases (at least away from the shock) . 

Shown in Figure 3.7 is the first-harmonic component of the unsteady shocked flow 
as computed from the time-marching code. The results are scaled by Ap ex it< Behind 
the shock, the unsteady pressure distributions are nearly identical for all of the cases 
although the case where Ap ex it is 0.080 differs very slightly from the other cues and the 
linear solution. Upstream of the shock, the flow is quiescent except for this last case. 
The shock passing upstream through the throat allows the downstream influence to be 
felt upstream of the throat. 

Table 3.1.2 shows the mean, first harmonic, and second harmonic components of the 
shock motion. The mew position of the shock is seen to move slightly forward with 
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Figure 3.6: Mean transonic flow in quasi-one-dimensional channel for various levels 
of unsteady back pressure. Note that the mean flow away from the shock is nearly 
indistinguishable from the steady flow in the channel. 
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Figure 3.7: First harmonic component of unsteady transonic flow in channel for various 
levels of unsteady back pressure. Note that the Fourier components for the different lev- 
els of back pressure are indistinguishable from one another except for the Ap eX jt = 0.08 
case. For this high level of unsteadiness, the shock passes through the throat, disap- 
pears, and then reforms aft of the throat. 
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increasing downstream back pressure variation. The situation is somewhat analogous 
to a nonlinear spring. It is easier to move the shock a finite distance forward than it 
is to move it aft. This occurs because the shock gets weaker as it moves forward and 
stronger as it moves aft. The first-harmonic content is seen to agree quite well with 
linear theory for the first two cases. The magnitude of shock motion agrees to within 
about 1.2 percent and the phase agrees to within less than 0.5 degrees. For the case 
of Apexit of 0.040, the error is about 3.5 percent in magnitude, while the phase is off 
by about 2.3 degrees. This is a truly remarkable result considering the total shock 
excursion is about 10 percent of the channel length for this case. For shock excursions 
much larger than this, the linear analysis is no longer valid as the shock disappears and 
reappears over the cycle. 

These results have important implications for the aeroelastic analysis of turboma- 
chinery. As a rough rule of thumb, a linearized analysis should be regarded as suspect 
if the amplitudes of the unsteady part of the flow flow is greater than ten percent of 
the total local values. For example, to predict the onset of flutter, which is governed 
by the behavior of the unsteady flow for very small flow perturbations, the unsteady 
flow may be modelled via a small perturbation analysis, other hand, the blades cannot, 
to predict the forced response of the blading due to incoming wakes or inlet distortion, 
the unsteady flow may or may not be analyzed accurately using a linearized analysis 
depending on the level of unsteadiness. 


3.2 Comparison to Full Potential Theory 

Some of the most widely used methods for predicting subsonic and transonic un- 
steady flows through real cascades are those based on linearized full potential the- 
ory [33,19]. However, the isentropic assumption inherent in the potential theory places 
a limit on its ability to accurately predict unsteady flows containing shocks (just as 
the inviscid assumption may place limitations on the linearized Euler method). In this 
section, we will compare results based on linearized Euler theory to those based on full 
potential theory. In particular, it will be shown that the isentropic assumption can 
produce serious errors in the predicted shock motion even for relatively weak shocks, 
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especially for low reduced frequencies. 

For some transonic flow problems, the full potential equations provide a useful ap- 
proximation. By assuming the flow is irrotational the velocity can be expressed as the 
gradient of a single scalar function, i.e., the velocity potential. If the flow is further 
assumed to be isentropic, the other primitive variables can also be expressed in terms 
of this potential. The equations which describe this irrotational and isentropic fluid 
motion are essentially a statement of mass conservation so that momentum and energy 
will, in general, not be conserved at shocks. Nevertheless, these equations have been 
shown to give reasonably good predictions of steady transonic flows for Mach numbers 
less than about 1.3. The purpose of this section is test whether linearized unsteady 
flow predictions based on the potential equations also give good results for such Mach 
numbers. 

Numerical Comparison of Full Potential and Euler Theories 

As a test of the linearized full potential theory, let us return to the example presented 
in Sections 2.2.7 and 2.3.7. We wish to compare the steady and unsteady transonic 
flows predicted by the Euler theories to those computed by the potential theory. First 
consider the steady flow problem. For the Euler solution, we specified the upstream total 
pressure and density and the downstream static pressure. Picking this back pressure 
fixes the location of the shock in the channel. Once the flow is choked, decreasing the 
back pressure further does not increase the mass flow through the channel. Further, the 
pressure at the exit of the channel is somewhat less than the pressure at the inlet due 
to the total pressure loss across the shock. The situation is somewhat different in the 
potential flow problem. If the flow is choked, then the pressure is a function of the local 
area alone. This function is multi-valued having a subsonic and a supersonic branch. 
But if the flow is choked, and the flow is subsonic at both the inlet and the exit, then (if 
the channel is symmetric) the pressure at the exit must be equal to the inlet pressure. 
In other words, one cannot arbitrarily set the back pressure and get a valid solution to 
the full potential equations. This is because there is no total pressure loss through the 
shock. Hence, there is no mechanism to set the shock position in the potential theory. 
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This situation is shown graphically in Figure 3.8. 

Therefore, to compare the potential and Euler solutions, we choose the back pressure 
of the potential solution to be consistent with the incoming static pressure, and the shock 
position to coincide with the Euler solution shock position. The Euler and potential 
base flows are shown in Figure 3.9. Note that ahead of the shock, the two flows are 
identical. Downstream of the shock, however, the pressure and density predicted from 
the potential theory are larger than those predicted from the Euler theory, while the 
velocity is smaller. If one were only interested in steady solutions, these differences 
might be considered small and acceptable. As we will see, however, these relatively 
small differences can produce much more significant errors in the predicted unsteady 
flows. 

Consider once again the unsteady transonic flow problem presented first considered 
in Section 2.3.7. Upstream, the total pressure and density are varied harmonically 
with a frequency u > of 1.257, while downstream the static pressure is held constant. 
Figure 3.10 shows the computed unsteady pressure distributions for the linearized 
Euler and linearized potential theories. Ahead of the shock, the two theories produce 
identical results. This is hardly surprising since ahead of the shock the flow is isentropic. 
Downstream, however, there are some differences between the two solutions. More 
disturbing, the computed shock motions differ significantly. The difference between the 
phases of the two predicted shock motions is 30.5 degrees. Furthermore, the shock 
displacement predicted by potential theory is 26 percent less that that predicted by 
Euler theory. Why then are there such large differences in the unsteady results? 

Part of the answer lies in the differences in the steady pressure, density, and velocity 
downstream of the shock. These quantities determine the speed at which pressure waves 
travel in the fluid. These pressure waves are sometimes referred to as the J~ and J + 
characteristics. The speeds at which these waves travel are given by 

cj- = U-C (3.25) 

c J+ = U + C (3.26) 

where C is the local speed of sound. Note that downstream of the shock, the character- 
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Figure 3.8: Effect of back pressure on the shock position of choked flow in a 

one-dimensional converging-diverging channel for Euler (upper) and potential (lower) 
theories. Note that for potential theory, the back pressure is unique, and cannot be 
varied to fix the shock location. 
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Figure 3.9: Steady base flow about which unsteady Euler and potential flows are cal- 
culated. Note the flows are identical ahead of the shock, but differ slightly behind the 
shock. 
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Figure 3.10: Unsteady transonic flow calculated from Euler and potential theories. 
Upstream the toted pressure is varied isentropically with a frequency w of 1.257 
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istic speeds differ for the Euler and potential theories. Hence, any disturbance which is 
reflected off of the exit boundary and propagates forward will move at the wrong speed 
in the potential theory. This upstream travelling pressure wave arrives at the shock at 
the wrong time, causing a phase error in the shock motion. This is at least part of the 
reason for the difference between the two theories. 

One could argue that this model problem is too contrived, that for more physical 
problems, a radiation boundary condition would be required at the exit flow boundary. 
Hence, this upstream travelling pressure wave would not be present (or at least less 
prevalent) and this problem would be alleviated. To some extent this is a fair criticism 
of this model problem. Therefore, for our next model problem, we will apply a radiation 
boundary condition at the exit. This is done by setting the incoming characteristic 
variable to zero. A characteristic analysis (see Chapter 5) shows that the incoming 
characteristic is given by 

for the Euler theory and 

wj- = u - j==pP = (1 ~ M)u + jw<f> (3.28) 

for the potential theory. The static-pressure exit-boundary conditions were replaced 
by requiring the incoming characteristic to be zero. Figure 3.11 shows the computed 
results. The agreement is somewhat better, but note that the predicted shock motions 
still do not agree that well. In particular, the difference shock motion predicted by the 
potential theory lags that predicted by the Euler theory by 26.8 degrees. This phase 
lag was identified as a problem with potential theory by Ashley [34]. 

Next, the two codes were run for the same case except that the forcing frequency 
was lowered to oj = 0.1. The results are shown in Figure 3.12 Note that now the 
disagreement becomes worse, especially for the predicted shock displacements. The 
magnitude of this shock displacement predicted using potential theory is 5.9 times larger 
than that predicted from the Euler theory. The phase difference between the predicted 
shock motions is 83 degrees. Earlier, we saw that the zero frequency shock displacement 
(i.e., the steady shock position) depends on the ratio of the back pressure to the total 
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Figure 3.11: Unsteady transonic flow with radiation outflow condition. Upstream the 
total pressure is varied isentropically with a frequency <jj of 1.257 
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Figure 3.12: Low frequency transonic channel flow (a; = 0.1). Note large difference 
between shock motions predicted by Euler and potential theories. 
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Amplitude of Shock Motion, \x t \ 



Figure 3.13: Frequency response of shock motion. Shown is the amplitude of shock mo- 
tion versus frequency of unsteady back pressure. Note that as the frequency approaches 
zero, the potential shock motion is unbounded while the Euler shock motion reaches a 
finite value. 

inlet pressure. In other words, the displ acement of the shock per unit change in this 
pressure ratio is equal to the slope of the steady shock position as a function of pressure 
ratio. For the Euler theory, this slope is finite. But in the potential theory, the steady 
shock can can occur anywhere behind the throat with the exit pressure fixed as shown 
in Figure 3.8. That is to say that the slope of the shock position curve is infinite . 
Hence, for low frequencies, the shock motion becomes very large and is unbounded as 
the frequency goes to zero. 

Figure 3.13 shows a frequency response plot showing the shock motion in the chan- 
nel versus the frequency of excitation. As before, the upstream total pressure and 
density are varied, while downstream, radiation boundary conditions are applied. Note 
that for the Euler theory, the shock motion amplitude approaches a finite value as the 
excitation frequency approaches zero. The potential theory, on the other hand, predicts 
an unbounded response at zero frequency. This phenomenon is due to the isentropic 
assumption of the potential theory which does not allow for a loss in total pressure 
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across the shock. Hence, there exists no mechanism to constrain the motion of the 
shock. It should be mentioned that it may be possible to modify the potential theory 
to partially account for the steady and unsteady entropy generation at shocks. Recent 
work by Fuglsang and Williams [35] demonstrated that such corrections greatly improve 
the steady and unsteady flow predictions of the transonic small potential methods. 

3.3 Summary 

In this chapter, some of the advantages and limitations of the present method have 
examined. It was demonstrated that the linearized unsteady Euler analysis compares 
well with the nonlinear Euler analysis so long as the unsteadiness in the flow is less 
than about 10 percent of the mean flow. Hence, the linearized Euler analysis is useful 
for most aeroelastic analyses. Furthermore, it is believed that the linearized Euler 
analysis provides a significant improvement over current linearized potential analyses in 
the modelling of flows with shocks. 
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Chapter 4 


Fundamentals of 
Two-Dimensional Linearized 
Euler Theory 


Unsteady flows present many problems in the design and operation of turbomachin- 
ery ranging from noise, to unsteady heat transfer, to stall and surge. In this work, 
however, we are mainly concerned with the aeroelastic aspects of unsteady flows. In 
other words, we wish to determine the unsteady loads felt by the engine structure due 
to the unsteady flow field. In particular, the two main aeroelastic phenomena of interest 
are the flutter and forced response of blade rows. For the purposes of studying these 
problems, we will assume that the flow is adiabatic, inviscid, and compressible, but 
may be rotational and nonisentropic. These assumptions are reasonable so long as the 
flow remains attached and the blade boundary layers and wakes are confined to narrow 
regions. 

In Section 4.1 the conservation principles for an adiabatic, inviscid, compressible 
fluid will be presented in both integral and differential form. These different forms of 
the fundamental physical laws which govern the flow are each useful for different parts 
of the unsteady flow analysis. For instance, the integral forms are preferred for the 
derivation of shock and wake jump conditions and for the numerical integration tech- 
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niques, while the differential form is useful in analyzing the far-field behavior of the 
fluid. In Section 4.2, the linearization of the nonlinear Euler equations is presented. 
These linearized equations are at the heart of the unsteady flow analysis. The lineariza- 
tion is carried about some nominal steady flow. Hence, before the unsteady flow can 
be analyzed, the steady solution must be known. In Section 4.3 the basic steady Euler 
solver is introduced. The steady Euler equations are solved using a Newton iteration 
procedure. Finally, in Section 4.4 the details of the linearized unsteady Euler solver are 
presented. This includes the discretization of the field equations, the linearization and 
discretization of the unsteady boundary conditions, and the numerical solution of the 
resulting equation set. 

4.1 The Conservation Laws 

The conservation equations jure well known (see for example [36]) but are discussed 
here for completeness. In this section, three different forms of the conservation equations 
are presented. Although all three are equivalent, their different forms are advantageous 
for different aspects of the analysis of unsteady flows. In Section 4.1.1 the integral 
form of the conservation laws will be presented for the case of a control volume fixed in 
space. In Section 4.1.2, this integral form of the conservation laws is extended to the 
case of a deformable control volume. This form is particularly useful for the derivation 
of shock and wake jump conditions. Finally in Section 4.1.3, the differential form of the 
conservation equations, known as the Euler equations are presented. These are useful 
for analyzing the far-field behavior of the unsteady flow. 

4.1.1 Integral Form of the Conservation Laws 

We will consider the conservation principles applied to a non-heat conducting, invis- 
cid fluid. The quantities to be conserved are mass, momentum, and energy. Consider a 
fluid passing through a region D enclosed by a surface S which is fixed in space as shown 
in Figure 4.1. The statement of conservation of mass is that the rate of increase of mass 
within the volume is equal to the rate at which mass flows into the region through the 
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Figure 4.1: Control volume used to derive momentum equations. The volume is assumed 
to be fixed in space with the fluid flowing through the volume’s surface. 

surface (the mass flux). In mathematical terms 

l o y t piD =-J 3 tV-ndS (4.1) 

The minus sign on the right hand side of Equation 4.1 is due to the fact that a velocity 
vector V in the direction of the outward normal n corresponds to mass flowing out of 
the region, not in. 

In a similar way, the conservation of momentum can be expressed. The statement 
of the conservation of momentum is that the rate of change of momentum of a group of 
fluid particles is equal to the sum of the body forces acting on the fluid particles plus 
the integrated surface forces acting on the surface. In the absence of body forces, we 
have that 

J D ^pVdD + f s pV{V-n)dS = - J s pndS (4.2) 

The first term in Equation 4.2 represents the rate of change of momentum within the 
region D , while the second term represents the rate at which momentum leaves the 
domain (the momentum flux). The two terms combined represent the rate of change 
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of momentum of the group of fluid particles which pews through the surface S at any 
instant of time. This then is equal to the sum of the forces acting on the boundary of 
the surface 5. 

Finally, the conservation of energy is considered. In this case, the net rate of energy 
increase of a group of fluid particles is equal to the work done on them by external 
forces. Hence we have that 

L + *r> dD + + 1r )v • = - /. » v < 4 - 3 > 

where t is the internal energy per unit volume of the fluid. The term pV 2 / 2 represents 
the kinetic energy per unit volume. The rate of change of total energy of a group of 
fluid particles is equal to the two terms on the left hand side of Equation 4.3. This 
in turn is equal to the rate of external work done by the static pressure acting on the 
surface 5. 

4.1.2 Integral Form of the Conservation Laws for Deforming Vol- 
umes 

The integral conservation Equations (Equations 4.1, 4.2, and 4.3) were presented 
for the case of a fixed volume in space. Although this form of the conservation laws 
is especially useful for deriving the differential form of the conservation equations, and 
consequently for the numerical integration of the conservation laws, a more general form 
of these equations allows for the domain D enclosed by the surface S to move rather 
than to be fixed in space. This latter form is useful for deriving the shock jump and 
wake jump conditions. 

Consider the situation shown in Figure 4.2. A surface 5 encloses a control volume 
D. The vector R traces out the surface of the control volume. This surface vector is 
a function of time, however, so that the enclosed volume is not constant. Hence, the 
flux through the boundary is not only a function of the fluid velocity V , but also of the 
velocity of the surface As before, the time rate of change of the enclosed volume 
of mass is equal to the mass flux through the surface, but the flux is now composed of 
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Figure 4.2: Nonstationary control volume. The flux through the control volume surface 
is composed of two parts: The flux due to the fluid motion, and the flux due to the 
motion of the control volume surface. 

two parts. 

Uo pdD= -IA y -li K ) nis 

Note that the time derivative appears outside the volume integral. 

Similarly, one can express the conservation of momentum and energy as 

s/»' Y " +//V(v-|*). ni s 

7t I D {e + pv2 M dD + / 5 ( e + py2 / 2 )( v “ ^ R ) nds = ~l s pV ' nds 

4.1.3 Differential Form of the Conservation Laws 

In the previous two sections, the conservations laws were presented in integral form. 
Often, however, it is more useful to express the conservation equations in differential 
form. The divergence theorem can be used to make this transformation. For example, 
consider the conservation of mass in integral form (Equation 4.1). Using the divergence 
theorem, the integral on the right hand side can be converted from a surface integral 


(4.4) 

(4.5) 

(4.6) 


80 



into a volume integral. After this transformation, the integral form of the conservation 
of mass becomes 

/ D + v ■ (pV)) JD = 0 (4.7) 

If Equation 4.7 is to hold for any arbitrary region of fluid, then it must be that the 
integrand vanishes identically everywhere within the fluid. 

The conservation of momentum and energy can be similarly reduced to partial dif- 
ferential equations. Together with the conservation of mass, these equations are the 
well known Euler equations which are given by 
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(4.8) 


where eo is the total internal energy per unit volume, and ho is the total enthalpy per 
uni t, mass. So at this point we have six unknowns and four equations. The unknowns 
are p, u, v, p, e, and ho. The additional equations needed come from the definition 
of enthalpy, and the relationship between the energy and the state of the fluid. These 
latter two relationships can be used to eliminate the total internal energy and total 
enthalpy from Equation 4.8 in favor of the four primitive variables p, u, v, and p. The 
enthalpy, and total energy are given by 


. _ («o + P) 

*° -~T~ 

«o = « + £P(« 2 + v J ) 


(4.9) 


(4.10) 


The internal energy e is in general a complex function of the state of the gas, but for an 
ideal gas with constant specific heats, the internal energy is given by e = ^ryP- Hence, 
under this assumption, the Euler equations can be expressed as 
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pv 

puv 

pv* + p 

^Lpv + \pv{ u 2 + v 2 ) 

With the Euler equations expressed in this form, we have four unknowns and four 
equations. This form, although slightly less general than Equation 4.8, will be used 
throughout this report. Implicit in these equations is the assumption of an ideal gas. 

Note that from Equation 4.8 we see that for steady flows we have that V* V/io is zero. 
This means that total enthalpy is constant along streamlines in steady flows. Hence, if 
the total enthalpy is constant in the upstream region of a given flow, the total enthalpy 
will be constant everywhere in the flow. In this case, the flow is said to be homoenergetic. 
If the flow is homoenergetic, the energy equation can be eliminated from Equations 4.8 
and 4.11, as well as one of the primitive variables. Eliminating the pressure p in favor 
of the remaining primitive variables, and making the same assumptions as before about 
the nature of the fluid gives 

pu pv 

, d 

pu 2 + p puv 

puv pv 2 + p 

whprp 

(«’+«’)] 

This form of the steady Euler equations simplifies the numerical computation of the 
steady flow. Having only three equations per computational node instead of four cuts 
the computational time by a factor of about two for the direct solver used in this 
research. 

4.2 Linearized Euler Equations 

In the first part of this chapter, the Euler equations were presented in a few of 
their many forms. In this section, the Euler equation are linearized to obtain the 
two-dimensional unsteady linearized Euler equations. As in Chapter 2, the analysis is 
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composed of two main parts. In the first part, the steady flow is determined. In the 
second part, the Euler equations are linearized about the steady operating point. The 
linearized Euler equations are then be used to compute the unsteady flow. As shown in 
Chapter 3, this unsteady solution will be valid so long as the resulting flow perturbations 
are relatively small. In this section, the linearized Euler equations will be developed. 

4.2.1 Field Equation Linearization 

As in Chapter 2, the full time dependent flow is then approximated as the sum of a 
mean and a small perturbation component. 


p(*, y,t) 

= ?(*»y) + ?(*»»><) 

(4.13) 

u{x,y,t) 

= tf(*,y) + «(*>y>0 

(4.14) 

v{x,y,t) 

= V(*»y) + v(**y>0 

(4.15) 

p(*»y> 0 

= ^(*,y) + p(*,y,0 

(4.16) 


Further, it is assumed that the perturbation is much smaller in magnitude than the 
corresponding mean component. It should be pointed out that this approximation is 
valid asymptotically as the perturbation in the flow tends toward zero. It is not an 
exact representation of the flow. As shown in Chapter 3, however, this approximation 
gives good predictions for moderate levels of unsteadiness. 

Next, the perturbation assumptions (Equations 4.13-4.16) are substituted into the 
Euler equations (Equation 4.11). Collecting terms of equal order and neglecting terms 
of higher order produces the equation for the steady flow, and the equations for the 
unsteady perturbation quantities. The steady flow equations are given by 

pu 1 [ py 

d pu* + p d puv 

— H 

dx puy dy py* + p 

PU + \pU(U 2 + V 2 ) ^PV + \pV{U 2 + V J ) 

which is just the original equations with time derivative terms set to zero. More inter- 


= 0 (4.17) 
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eating are the first order perturbation equations, which are given by 

1 0 0 0 

U p 0 0 

V 0 p 0 

+ pU pV ^r 

U p 0 0 

U 2 2 pU 0 1 

UV pV pU 0 

\ U(U 2 + V 2 ) JL-P + *pU 2 + IpV* pUV £ 

VO p o' 

d_ UV pV pU 0 

d y V 2 0 2 pV 1 

\V(U 2 + V 2 ) pUV ^P + \pU 2 + \pV 2 $=i 
These equations are linear in the perturbation variables. The mean flow components 
are known from the solution of Equation 4.17. Hence, the matrices premultiplying 
the perturbation variables are known variable coefficient matrices. Note also that the 
perturbation equations are in so-called conservation form. One could in principle time 
march these linear equations to obtain an unsteady solution to the Euler equations valid 
for small deviations away from the mean flow. However, if the flow of interest is periodic 
in time, then the unsteady perturbation can be represented as a Fourier series. Hence, 
for example, the density can be expressed as 

+00 

p{x,y,t)= £ p n {x,y)e inut (4.19) 

n=-oo 

where u; is the fundamental frequency given by w = 2 x/T. But since the first-order equa- 
tions are linear, the behavior of each Fourier component can be analyzed individually, 
then summed together to form the total solution. Hence, without loss of generality, the 
Fourier components will be analyzed term by term. Furthermore, many flows of interest 
are not only periodic, but are also nearly harmonic. Hence, we let 

p(x,y,t) -* p(x,y)e jut (4.20) 
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u(*,y,t) -» u(x,y)e iui (421) 

v(*,y,0 - v(x,y)^* (4.22) 

p{x,y,t) -* p(*» y)e^"* ( 4 - 23 ) 

The perturbation amplitudes p, u, v, and p are complex. Substitution of the harmonic 
assumption into the linearized Euler equations gives 
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(4.24) 


[ \V{U 2 + V 2 ) pUV P + \pU 2 + | pV 2 
Equation 4.24 can be expressed in matrix notation as 

d d 

jwBiU + — B 2 U + — B 3 U = 0 (4-25) 

ox dy 

Thus we can solve for all time by considering a single set of linear variable-coefficient 
equations for each frequency of interest. Since only a single set of equations must be 
solved, the linearized approach should be computationally more efficient than one in- 
volving time marching. The former is well-suited for the flutter problem in cascades as 
well as the inlet distortion problem. It is not so well-suited for the case of an unsteady 
flow through a cascade which is driven by the wakes from the upstream stator be- 
cause this sort of disturbance, although periodic, is not harmonic. This problem could, 
perhaps, be better handled by time marching. However, one could in principle decom- 
pose the wake disturbance into Fourier components, solve for each of these components 
individually, and then sum the results to obtain the total solution. 
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4.3 Solution of the Steady Euler Equations 


In the previous section, the linearized unsteady Euler equations were presented. 
Because these equations are linearized about some nominal steady flow, the steady 
flow must first be calculated before the unsteady analysis can proceed. In this section, 
the basic method for calculating the steady flow is discussed. This procedure closely 
parallels the procedure used to calculate steady one-dimensional flows first presented in 
Chapter 2. The nonlinear steady Euler equations are solved using a Newton iteration 
procedure. The nonlinear problem is reduced to a series of linear partial differential 
equations, each very similar to the unsteady linearized Euler equations. Each of these 
linear problems is discretized on a computational grid and solved directly. The details 
of this entire process are given below. 

4.3.1 The Linearized Newton Iteration Equations 

The Newton iteration method for solving two-dimensional partial differential equa- 
tions is an extremely fast algorithm [24]. Furthermore, the method reduces the non- 
linear steady Euler equations to a series of linear ones, each of which is similar to the 
linearized unsteady Euler equations. Each linear problem can be solved directly using, 
for example, a Gaussian elimination routine. Unfortunately, the Newton method is not 
guaranteed to converge. With some care, however, the procedure usually converges. 

Consider the steady Euler equations (Equation 4.17). Let us assume that an estimate 
of the actual solution is known. Then the exact solution is assumed to be the sum of 


the approximate solution plus a small correction. Hence 


p(x,y) = p(x,y) + p(x,y) 

(4.26) 

u{x,y) = U{x,y) + u(z,y) 

(4.27) 

v(x,y) = V(x,y) + v(x,y) 

(4.28) 

p(x,y) = P[x,y) + p(x,y) 

(4.29) 


where the left hand side is the “exact” solution, or at least a much improved estimate. 
The two terms on the right hand side are the current estimate of the solution and the 
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correction, respectively. The correction is assumed to be small compared to the current 
estimate of the solution. Substitution of the assumed solution into the steady Euler 
equations, expanding the result into a power series, and then truncating second- and 
higher-terms gives the steady linearized Euler equations, i.e., 


d „ d _ (dY , 3G\ 

— b 2 u + 3 y BsU Ux + a»j 


(4.30) 


where u is the vector of correction perturbations The term on the right-hand 

side are the residuals of the steady equation based on the current estimate of the flow 
[p,U,V] T . The left hand side represents the first order correction to the residuals due 
to the correction perturbations u. The variable coefficient matrices are given by 


B, = £ (4.31) 

B, = §§ (4.32) 

Note the similarity of the Newton iterations equations to the linearized unsteady Euler 
equations (Equation 4.24). The two linearizations are basically the same. The main 
differences are that the linearized unsteady equations contain a homogeneous term aris- 
ing from the time-derivative term, while the Newton iteration equation contains an 
inhomogeneous term which represents the residuals of the current solution estimate. 

The basic Newton iteration procedure is to start with a reasonable guess for the 
steady flow solution. The iteration equations are then discretized on a computational 
grid and solved subject to an appropriate set of boundary conditions. Having determined 
the corrections, these are then added to the current estimate of the solution to obtain an 
improved estimate. The entire process is repeated until a converged solution is obtained. 
For subsonic flows, convergence usually occurs after about five iterations. For transonic 
flow calculations, the procedure usually converges in from five to ten iterations. 


4.3.2 Discretization of the Newton Iteration Equations 

The Newton iteration equations are linear variable coefficient equations. In gen- 
eral, however, these equations cannot be solved exactly and, therefore, must be solved 
numerically. In this section, the discretization procedure is discussed. 
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Figure 4.3: Typical grid used in calculation of steady and unsteady channel flows. 


88 


A typical grid used to calculate channel flows is shown in Figure 4.3. The nodes 
are referenced by two subscripts («, j) denoting the node number in the streamwise and 
normal directions. In the streamwise direction, the nodes are numbered from 1 to /, 
while in the normal direction the nodes are numbered form 1 to J. There are a total of 
(/ — 1)( J - 1) cells. At each of these cells, the Newton equations are discretized using a 
finite volume operator similar to the one used by Ni [23], Consider once again the New- 
ton equations. If the Newton equations are satisfied everywhere within a conservation 
cell, then the area integral of the equations is also satisfied and given by 

// (£ b ’ u + £■*») did >=- // (I + w ) d * d » <« 3 > 

where the limits of integration are over a single conservation cell. Using the divergence 
theorem, this area integral is converted into a line integral of the form. 

£ (Bju.Bati) • n da = - ^(F,G) • n da (4.34) 


where n is the outward unit normal of the cell. This is the form of the Newton equations 
which are discretized. 

Consider a typical conservation cell as shown in Figure 4.4 The primitive variables 
are stored at the cell corners. That is to say the scheme is a node-centered scheme. The 
trapezoidal rule is used to integrate the line integrals in Equation 4.34. For example, to 
calculate the right hand side, the values of F and G are calculated at each of the nodes. 
These values are averaged from the two nodes on each face, the dot product with nAs 
is formed, and the contributions from the four faces are summed. Hence, integration 
around the cell gives 


/(F,G) • nda » + 5f±Lii^Wti Ayj 

+ r. t u t i + ?wi a|[ + f m Av4 _ ?» + < ■> +& 

_ G <+i„> + G «- +1,;+1 Ax2 _ G «-n,i+i + G «.y+1 Ax& _ G «•.;+ 1 + G «'.i Ax4 

2 2 2 


(4.35) 


where 


Al! = Xj + i,jf - Xij 
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Figure 4.4: Discretization of residuals using a finite volume integration 
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AX2 


*i+X,J+X ~ x i+l,j 

Azj 

= 

x i,i+l ~ *i+l.i+l 

Ax 4 

= 

x i,j ~ x i,j+l 

Ayx 

= 

y»+i.i ~ vi,i 

Ay2 

= 

y.+x,y+i - y»+i,y 

Ays 

= 

Vi,i+l - Vi+l.j+l 

Ay 4 

= 

vi,j - y«-,j+x 


Rearranging slightly gives 


Ayx + Ayj 


f (F, G) • n da w F t> , Ay4 + Vl + F t+1> ,- 

. „ Ayj + Ays , „ Ays + Ay 4 „ Ax 4 + Az x 

+* «+x,i+i j + * i,i+l 2 '* <J 2 

„ An + A®2 ^ A*j + A®s „ Azj + Ax 4 

- G.+ i,i j 2 2 

Similarly, the homogeneous portion of Equation 4.34 is discretized by 


(4.36) 


<J> (B 2 u, B311) • nds w (B 2 u),,j Ay4 y Ayi - + (Ban)<+ij — y — — 


x Ayj + Ays . 4 Ay s + Ay 4 m ^ A z 4 + An 

+(Bju),+x iJ +i h (B2 u), iJ+ x ^ (Bsu),j ^ 

„ An + Az 2 x Az 2 + Ax s ra ..x Az 8 + Az 4 u ^ 

- (Bsujj+x,,- (B s u)i + x t y + x (B 3 u)i J+1 ( 4 . 37 ) 


This discretization scheme is second-order accurate. Furthermore, because the scheme 
is node centered, no extrapolation to the boundaries is required, as with cell-centered 
schemes, to apply boundary conditions. The major flaw with the scheme is that it 
admits sawtooth or hourglass modes. A sawtooth mode averages to zero at the face 
centers. Hence, such modes do not contribute to the numerical integration around the 
cell. To eliminate such modes, a small amount of smoothing must be added to the 
scheme. This is true of Ni’s time-marching scheme as well. The addition of smoothing 
is not totally undesirable as some smoothing must be added to stabilize the solution at 
the sonic line for transonic problems. Regardless, very small amounts of smoothing are 
required in practice, but not enough to significantly alter the solution. 
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4.3.3 Steady Flow Boundary Conditions 

In addition to the discretization of the conservation equations, appropriate boundary 
conditions must be specified, linearized, and then discretized for use in the Newton 
iteration procedure. To start, consider the boundary conditions associated with two- 
dimensional channel flows. For channel flows, there are three boundary regions: the 
inflow boundary, the outflow boundary, and the channel walls. At each of these regions, 
one or more boundary conditions will usually be applied. 


Flow Tangency Condition 

Along the solid wall boundaries of the channel, the flow is required to be tangent to the 
wall. That is to say that no mass flows through the wall. This boundary condition is 
expressed as 

(u,t5) • n = 0 (4.38) 

where n is the local unit normal to the wall. Expanding Equation 4.38 in terms of the 
current estimate of the flow and the correction, we have that 


(u,v) • n = -(17, V) • n 


(4.39) 


As this flow tangency condition is already linear in the correction perturbation variables, 
no linearization step is required. 

This boundary condition is applied at the face centers to the cell faces along the 
channel walls. The normal is taken to be the normal of the face center and the velocities 
are averaged from the nodes at the ends of the face. Hence, the discretized flow tangency 
condition is given by 


where 


A* — x t+i,j x i,jt Ay — yi+ij yi,j 
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Outflow Boundary Conditions 

For the channel flows considered in this report, the outflow is subsonic. Hence, a 
single boundary condition is applied along the outflow boundary. Although a variety 
of boundary conditions could be applied, in the present effort, the static pressure will 
be specified at the outflow boundary. Since the steady Euler equations are being solved 
with a reduced set of primitive variables, the pressure must be specified in terms of the 
density and velocity. Hence, the exit pressure is given by 

Petit = 1 ^ 1 [pho- \p («’ + t> 1 )] (4.41) 

Linearizing this boundary condition about the current estimate of the flow gives 

[to - j (ff* + V>)] , - = 

(17* + V*)] 

= Petit - P (4.42) 

Equation 4.42 is applied at each of the J nodes along the outflow boundary. 


Infl ow Boundary Condition 


For a subsonic homenergetic inflow, two boundary conditions are required at the inflow 
boundary. In the present work, the flow angle and total density will be specified. A third 
implicit assumption due to the homoenergetic assumption is that the total pressure is 
specified. The total density is given by 


PT in let-p\l 2hQ J 


which, when linearized, becomes. 


(7 - 1)^0 


U* + V 2 ' 

2 ho 


PV_ 

(7 " l)*o 


u* + v 2 

2/10 


j v = PT inlet - P ^ “ 


2ho 


Equation 4.44 is applied at each of the J nodes along the inflow boundary. 


93 


The last boundary condition to be specified is the angle of the incoming flow. For a 
channel whose walls are parallel at the inlet, the flow at the inlet will be parallel to the 
walls. For channel whose walls are not parallel, the specified flow angle will vary along 
the inflow boundary. Nevertheless, at each node on the inflow boundary, the flow angle 
must be specified. Let this flow angle be denoted by a. Then 

usin(a) - t)cos(a) = 0 (4.45) 

or 

usin(a) - vcos(a) = -U sin(a) + V cos(a) (4.46) 

Equation 4.46 is specified at each of the inflow boundary nodes. 

At this point, we should determine whether the number of equations is enough to 
specify all of the unknowns. Figure 4.5 shows where the conservation equations and 
boundary conditions are applied on a grid used to calculate flow in a channel. Since the 
computation grid has IJ nodes, the total number of unknowns is 37J. The equations 
available to solve for these unknowns are the conservation equations plus the various 
boundary conditions. There are (/ — 1)(J — 1) conservation cells with 3 conservation 
equations per cell for a total of 3(7 — 1) (7 - 1) equations. There are 2 J inflow boundary 
conditions, J outflow boundary conditions, and 2(7— 1) wall boundary conditions. The 
total number of equations is then 37J — 7+1 equations. In other words, there Me 7—1 
equations too few. This shortage of equations is due to the method of discretizing the 
equations. These 7—1 degrees of freedom are sawtooth modes admitted as solutions to 
the rank deficient system of equations. To eliminate them, smoothing will be added to 
the equations after preconditioning them as discussed in the next section. 

4.3.4 Matrix Equation Structure 

Written in shorthand form, the discretized equations which describe the linearized 
unsteady flow can be expressed as 

Mu = f (4.47) 

The matrix M is a supermatrix which is block bidiagonal. Each of these blocks is in 
turn also a supermatrix which is block bidiagonal. Each of these blocks is a three by 
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Figure 4.5: Location of conservation equations and boundary conditions on a typical 
channel flow computational grid. 
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three matrix. Note in particular that the last row of the supermatrices is zero. This is 
because there are IJ nodes but only (/ — 1 )(J — 1) conservation equations. 

Similarly, the boundary conditions can be expressed in matrix form as 

Mbcu = f bc (4-48) 

where the matrix M bc has many fewer rows than columns and is sparse. 

An equivalent system of equations can be formed by "squaring” the conservation 
equations and the boundary conditions and summing them to give 

[M r M + M£ c Mj,c] u = M r f + *A T BC t BC (4.49) 

One interpretation of Equation 4.49 is that the solution of it minimizes the square of the 
residuals (sets to zero) of the original equation set. Of course, the “squared” equations 
are still rank deficient and admit sawtooth modes. However, the matrix equations are 
now positive semi-definite. The rank deficiency is eliminated, along with the sawtooth 
modes, by adding a small amount of positive semi-definite smoothing to the equation 
set. A finite element Laplacian operator L is used. Hence, the resulting system of 
positive definite equations is given by 

[M T M + M£ c Mbc + «L] u = M r f + Ml c i bc + di (4.50) 

where 

f L = -LU 

and c is a small parameter introduced so that only a small amount of smoothing is added 
so as to eliminate the sawtooth modes without significantly modifying the solution. 

This process of forming the equations to be solved can be viewed as a kind of 
constrained optimization. The quantity to be optimized is the smoothness of the solution 
subject to the constraint that the conservation equations and the boundary conditions 
be satisfied. This construction is similar to techniques in numerical optimizations. A 
cost function to be minimized is modified by adding constraint penalty functions. For 
this example, the modified cost function is given by 

R = i (Mu - f) 3 + i (M B c« - tBcf + J(u + U) T L(u + U) (4.51) 

2 2 2 
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Minimizing this scalar with respect to the perturbation variables u gives 

^ + cLj u — (M r f + M % c*bc + eft) = 0 (4-52) 

which is exactly equivalent to Equation 4.49. 

This provides a rational approach for the addition of smoothing to the problem. 
Furthermore, the resulting equations are positive definite. Unfortunately, the resulting 
equations tend to be poorly conditioned. Hence, double precision arithmetic is used 
in the computer code to avoid potential round-off problems when performing Gaussian 
elimination on the matrix equations. The bandwidth of the matrix is also increased. The 
super matrices are now block tridiagonal. Said another way, the original equations were 
derived from operations which involved four nodes. The new operator is a nine-point 
operator. 

The foregoing matrix equations are solved directly using a Gaussian elimination 
scheme. The solver takes advantage of the known sparseness, but does not reorder the 
equations for maximum efficiency. The computer time required to solve this system of 
equations is 0(/(3J) 8 ). The amount of storage required is 0(/(37) J ). Hence, addi- 
tional grid resolution in the normal direction greatly increases both CPU and storage 
requirements. 


4.4 Unsteady Flow Calculations 

In the previous section, the numerical method for solving the steady two-dimensioned 
flow through a channel was described. The steady flow must be calculated before one 
can solve the unsteady Euler equations to obtain unsteady flow predictions. In this 
section, the numerical method for solving the linearized unsteady Euler equations will 
be detailed. In Section 4.4.1, the method of discretizing the linearized unsteady Euler 
equations is discussed. In Section 4.4.2, the treatment of the boundary conditions for 
unsteady channel flow are presented. The three steps involved are the specification of 
the (nonlinear) boundary conditions, the linearization of the boundary conditions, and 
the discretization of the linearized boundary conditions. Finally, in Section 4.4.3, the 
conditioning of the matrix equations is discussed. 


97 


4.4.1 Linearized Unsteady Euler Equation Discretization 


The discretization of the linearized unsteady Euler equations is similar to the dis- 
cretization of the Newton iteration equations. The starting point is the linearized un- 
steady Euler equations, i.e., 

Q Q 

jwBiU + — B 2 u + — B s u = 0 (4.53) 

If these equations hold over an entire conservation cell, then the integral over the cell 
must also be zero. Hence 


J I (ywBiU + ^B 2 u + j^Bju) dx dy = 0 (4.54) 

The second and third terms can be converted to a line integral by application of the 
divergence theorem. The resulting integral equation is given by 


J J ju Biu dxdy + j (B 2 u, B 2 u) • n da = 0 


(4.55) 


The next step is to numerically aproximate these integrals. The last two integrals can be 
evaluated using a trapezoidal integration around the cell as with the steady equations. 
The first term, however, is approximated by the average value of the integrand at 
the four corner nodes multiplied by the area of the conservation cell. The discretized 
linearized Euler equations are then given by 


J j jwBiU dxdy + j (B 2 u, B s u) • n da fa 


+(Bm)w 


(Biu)i,j + (Biu),+i,j + (Biu), •+!,,+! + (Bxu), >J+ i 

4 

Ayi + Ajfe 


juiA 
Ay 4 + Ayi 


+ (B 2 u) i+ i ; - 


+ (B 2 u),+i >j+ x 


Ajfc -I- Ays 


- (B,u) i+1 , /+1 t — - t = 0 < 4 ' 56 ) 


where A is the area of the ijth cell. This discretization of the linearized Euler equations 
is second-order accurate in the linear dimension of the cell. Note further that the 
equations are homogeneous and the perturbation variables are complex. 
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As with the discretization of the steady Euler equations, the discretization of the 
linearized unsteady Euler equations admits sawtooth modes. This problem is typical of 
node based schemes. For this reason, a small amount of smoothing must be added to 
the assembled equation set. 

4.4.2 Boundary Conditions of the Unsteady Flow Solution 

XJp to now, nothing has been said about the boundary conditions which will be 
required to specify the unsteady flow in a channel. These conditions may reflect a real 
physical problem to be solved, or may be somewhat more contrived. Regardless, they 
must be well posed. For channel flows, the boundary conditions can be divided into 
three regions of the flow: inflow, outflow, and solid wall boundary conditions. The 
number of boundary conditions required at the inflow and outflow boundaries depends 
on the local Mach number. As with the linearized unsteady Euler equations, setting 
up the numerical equations is divided into three parts: specification of the governing 
boundary conditions, linearization, and discretization. 

Solid Wall Boundary Conditions 

Consider the wall boundary shown in Figure 4.6. Because the wall is solid, there is 
no mass flux through the surface. This is the flow tangency condition as expressed 
mathematically as 

(u, v) • n = 0 (4.57) 

where n is the outward normal to the wall surface. Expansion of Equation 4.57 in terms 
of the mean flow and the unsteady flow components gives 

(ue* w ‘, • n = -{U, V) • n (4.58) 

But because the mean flow satisfies the flow tangency condition, the right hand side of 
Equation 4.58 is zero. Hence, the unsteady flow tangency condition is given by 

(u,v)-n = 0 (4.59) 

Because this boundary condition is already linear, the linearization step can be skipped. 
The fined step is the discretization of Equation 4.59. The boundary condition is applied 
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at the centers of the wall boundary faces. The velocity at the face center is taken to be 
the average of the velocities at the two nodes, and the normal is taken to be the normal 
of the face. Hence, the flow tangency condition is discretized as 

_ u ».i + Ay + V JA + Vi ±hi A x = 0 (4.60) 

2 2 


where 


Az = *t+ i,j — *«,;» Ay — y«+i,y — y«,/ 


Outflow Boundary Conditions 

For the two-dimensional channel flows considered in this report, the flow at the exit of 
the channel is subsonic. Hence, one boundary condition is required. For the moment, 
the unsteady back pressure is specified along the outflow boundary. The unsteady back 
pressure at the exit of the channel is given by 

p = Petit + Ap«*«e ,w * (4.61) 

where P ex n is the steady exit pressure. Hence, the unsteady portion of the outflow 
boundary condition is given by 

P = Ap e «tf ( 4 - 62 ) 

This boundary condition, like the flow tangency condition, is also linear so that the lin- 
earization phase of the discretization of the boundary condition is not required. Equa- 
tion 4.62 is applied at each of the nodes along the outflow boundary. 

Note that this is a reflecting boundary condition. That is to say that waves which 
impinge on the outflow boundary will be partially reflected. Thus Equation 4.62 is a 
nonphysical boundary condition. In the next chapter, a more sophisticated nonreflecting 
boundary condition will be developed for cascade flows. 

Inflow Boundary Conditions 

For subsonic inflow, three boundary conditions must be specified at the inflow boundary. 
For now, the quantities specified will be the total pressure, total density, and inlet flow 
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angle. The unsteady total pressure variation is assumed to be 


Pr = Pro + A Pre*"* 


(4.63) 


where Pro is the mean total pressure and A Pr is a small perturbation in the total 
pressure. Similarly, the total density is given by 


Pt = pro + A prt* ut 


Hence, we have that the total enthalpy is given by 
1 Pr 7 Pro + A Pre* wt 


7~1 Pr 7 — 1 Pro + A pre iwt 7 

Similarly, a quantity closely related to the entropy is 

Pr p Pro + A Prc’ ut 




Pt 


(4.64) 


(4.65) 


(4.66) 


P (pro + Apre Jwt y 
These boundary conditions are clearly nonlinear functions of the total pressure and 
density perturbations. Expanding the expressions into perturbation series in powers of 
APr and Apr and truncating the second- and higher-order terms gives the zeroth- and 
first-order boundary conditions. The zeroth-order equations are satisfied by the mew 
flow. The first-order unsteady boundary conditions are 

A Pr 7 Pro 


"^T G) p - (£) ,+t '“ +Vv= T^T 


Pro 7 - 1 Pro 


and 


P 

p p_ 

Pro APr Apr 

p 

p 

Pro -Pro Pro . 


Apr (4.67) 


(4.68) 


Note in particular that if the flow is isentropic to first-order, then the right hand side 
of Equation 4.68 (and hence, the left) is zero. 

The inlet flow angle is assumed to be time-invariant for the channel flow problems 
considered in this report. Hence 


usin(a) - t)cos(a) = 0 


(4.69) 


which when linearized becomes 


u sin(a) - v cos(a) = 0 


(4.70) 
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4.4.3 Matrix Equations 

The assembly of the linear equations for the unsteady analysis is similar to that for 
the Newton equations used in the steady flow analysis. Because the unsteady equations 
are complex, however, the transpose operator is replaced by the Hermitian operator. 
The resulting matrix equations are 

[M*M + Mg c M BC + cL] u = Mg c f BC (4.71) 

Because the linearized unsteady Euler equations are homogeneous, the only inhomoge- 
neous term in Equation 4.71 comes from the boundary conditions. This matrix equation 
is Hermitian and positive definite. 

4.5 Summary 

In this chapter, the basic method used to numerically calculate two-dimensional 
unsteady flows has been presented. The major steps in the analysis are the discretization 
of the governing nonlinear equations, the linearization of the boundary conditions, the 
discretization of the linearized field equations and boundary conditions, the assembly of 
the discretized equations into matrix form, and the solution of the matrix equations. At 
this point, fairly simple two-dimensional flows may be analyzed. In the next chapter, 
the extensions of this method which are needed to analyze more complicated flows, such 
as cascade flows or flows with shocks, will be presented. 
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Chapter 5 


Extensions of Linearized Euler 
Theory to Cascade Flow 


In Chapter 4, the fundamental aspects of the linearized Euler analysis were pre- 
sented. With this foundation, fairly simple two-dimensional channel flows can be ana- 
lyzed. In this chapter, the extensions to the theory which are needed to analyze more 
complicated steady and unsteady cascade flows are presented. In particular, the exten- 
sions needed to accommodate steady and unsteady shocks, steady and unsteady wakes, 
moving airfoil boundaries, and acoustic radiation in the far-field are addressed. The 
treatment of these flow features closely parallels the methods used by Verdon [33] and 
Verdon, Adamczyk, and Caspar [37] to analyze potential flows. 

Shown in Figure 5.1 are the five main boundary surfaces to be considered in this 
chapter. They are the moving airfoils, the upstream periodic boundaries, the down- 
stream wake slip-planes, the upstream far-field boundary, and the downstream far-field 
boundary. Furthermore, if a shock occurs in the domain, a sixth internal boundary must 
be considered to connect the two regions of continuous flow upstream and downstream 
of the shock. In Section 5.1 the boundary conditions applied at moving airfoils will be 
considered. This boundary condition will allow the flutter problems to be analyzed. 
In Section 5.2 the treatment of steady and unsteady shocks and wakes is discussed. 
As in the previous one-dimensional analysis, a fitting procedure is used to accurately 
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describe the mean position and the motion of the shock. Similarly, wake fitting is 
used to model the mean position and motion of the wake. In Section 5.3, the periodic 
and far-field boundary conditions will be considered. The use of periodic allows the 
steady and linearized unsteady flows to be analyzed by looking at a single blade pas- 
sage, greatly reducing the required computational effort. Finally, the treatment of the 
far-field boundary conditions is examined. Because the computational domain must be 
finite, the far-field boundaries are placed at a finite distance in front of and behind the 
cascade. Yet they must allow disturbances to pass through them as if they were not 
there. Hence, special nonreflecting boundary conditions must be employed. 

5.1 Moving Airfoil Boundaries 

In this section the boundary conditions applied to a moving airfoil are discussed. 
Although the analytical boundary condition is straightforward, the numerical boundary 
condition is complicated slightly by the fact that the grid is fixed in space while the 
airfoil is supposed to vibrate harmonically. 

To start, consider the boundary condition for a solid airfoil moving through a fluid. 
Suppose the surface of the airfoil at any point in time is described by the position vector 
R. Then the boundary condition (flow tangency) at the surface of the airfoil is 

V • n - • n = 0 (5.1) 

where n is the unit normal to the surface. This boundary condition simply states that 
there is no mass flux through the surface of the airfoil. Note further that the unit normal 
to the airfoil n is a function of the airfoil surface position vector R. Hence, Equation 5.1 
is nonlinear. To linearize this boundary condition we assume that the airfoil surface 
position, like the field variables, can be represented as the sum of two components: one 
describing the mean location and the other a small perturbation in this location which 
is a function of time, i.e., 


R(s,t) 

= R(s) + r(s,t) 

(5.2) 

n(s,t) 

= ff(«) + n(«,t) 

(5.3) 
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f(a,t) = r(«) + r(a, t) (5-4) 

where t is the unit vector tangent to the airfoil surface and the airfoil surface vector R 
has been expressedas a function of the distance s along the surface of the airfoil. Then 
the unit normal n and unit tangent t at the moving airfoil surface are related to the 
corresponding vectors B and f and the mean surface position by 

h = B-r^ (5.5) 

08 

f = r + K 5=- (5.6) 

da 

These equations are asymptotically valid for small airfoil motions. 

Substitution of the perturbation assumption for the motion of the airfoil surface 
and the primitive flow variables into the flow tangency condition (Equation 5.1), and 
collection of terms of equal order gives the zeroth-order (mean) and the first-order 
(linearized unsteady) boundary conditions. The mean-flow boundary conditions are 
simply 

V • H = 0 (5.7) 

The Newton linearization and discretization for use with the steady solver was discussed 
in Chapter 4. 

The first-order unsteady boundary condition is 

'■*=ir +v ‘ a i? (58) 

where the subscripts n and t denote the component in the normal and tangential direc- 
tions respectively. Equation 5.8 says that the upwash on the blades is composed of two 
components. The first is due simply to the velocity of the airfoil normal to its surface. 
The second term is due to the rotation of the airfoil. If the flow is quasi-steady so that 
the first term disappears, an upwash will still be induced by rotation of the airfoil. A 
nose down motion of the airfoil requires that a positive upwash occur so that the flow 
will remain tangent to the airfoil surface. 

It should be emphasized that Equation 5.8 must be applied at the instantaneous 
position of the airfoil surface and not the mean location. This is a problem computa- 
tionally since the grid we use is fixed in space. Hence, an additional term is required to 
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extrapolate this boundary condition from the instantaneous airfoil position to its mean 
position. As a result, the boundary condition at the mean surface is 

V 'H= ^ + - r - v (V„) (5.9) 

If the flow to be analyzed is harmonic, then the operator is replaced by ju>. 

Equation 5.9 can be discretized in a straightforward manner. The numerical bound- 
ary condition is applied at the centers of the cell faces which lie along the boundary 
of the airfoil surface. The mean normal and tangent of the airfoil surface are taken to 
be the normal and tangent to the cell face. Because the airfoil motion is prescribed, 
r, r n and are known. The only real difficulty is in the evaluation of the gradient 
of the normal component of the mean velocity V„. This is evaluated using a first-order 
accurate difference operator which uses the values of the mean normal velocity V n at 
the four corners of the boundary cell. Because this procedure errors of O(Ax) into the 
scheme at the airfoil boundaries, there is the possibility of degrading the accuracy of 
the entire method. Fortunately, the cells next to the wall are usually very small, and 
the gradient term is often small, limiting the effect of this first order error. Although 
not implemented at present, a second order approximation could be used. 

5.2 Flow Discontinuities 

The Euler equations are nonlinear and admit so-called weak solutions. Weak so- 
lutions are those which are not everywhere differentiable, but nonetheless satisfy the 
integral form of the conservation equations. Flows which contain shocks or wakes fall 
into this category. In real flows, shocks and wakes have some small but finite thickness. 
However in the absence of viscosity, shocks and wakes are modelled as surfaces at which 
the flow variables are discontinuous. In this section, the so-called jump conditions which 
govern the behavior the flow at these surfaces are developed. These jump conditions 
will subsequently be used to fit shocks and wakes in the present steady and unsteady 
flow analyses. 

In Sections 5.2.1 and 5.2.2, the steady and unsteady wake- and shock-jump condi- 
tions will be discussed. These conditions relate the flow on one side of a wake or shock 
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Figure 5.2: Control volume used for the derivation of the shock jump conditions. The 
volume goes to zero in such a way that only a small portion of the shock is enclosed. 

to the flow on the other side. These nonlinear jump conditions will be linearized for for 
use with the linearized unsteady analysis. 

5.2.1 Unsteady Shock Jump Conditions 
Nonlinear Shock Jump conditions 

To derive the unsteady shock jump conditions, we make use of the integral conservation 
laws as expressed Equations 4.4, 4.5, and 4.6. Consider a volume D which contains a 
segment of a shock as shown in Figure 5.2. The shock position is described by the vector 
R. The vectors n and r are the unit normal and the unit tangent to the shock surface, 
respectively. In general, the shock will not be stationary, but will move through the 
fluid. Next, we imagine the surface 5 to shrink until it just encloses a small length of the 
shock. Evaluations of the integrals in Equations 4.4, 4.5, and 4.6 gives the desired shock 
jump conditions from the integral conservation equations. For example, the continuity 
of mass requires that 


dR 



The symbol [[• • •]] denotes the difference in the enclosed quantity across a discontinuity. 
Equation 5.10 enables us to connect the two regions of continuous flow on either side of 
a shock. Similarly, the momentum and energy conditions are given by 

l(pVV • n + pn]| — f[pV|| • n (5.11) 

[[(« + p + ^ ! /2)V.n]]= [[(« + ^V2)]]||-n (5.12) 

Note in particular the important role the normal component of the shock velocity has in 
the shock jump conditions (Equations 5.10-5.12. In fact, the tangential direction does 
not appear at all. Without much loss in generality, we can think of the motion of the 
shock as being normal to its surface. Alternatively, we can consider the shock to move 
in the direction of some arbitrary line (so long as this line is not tangent to the shock) 
such as a computational grid line since there is no constraint on the tangential motion 
of the shock. 

The momentum shock jump condition can be further simplified by looking at the 
components in the direction normal to and tangential to the shock. The normal and 
tangential momentum shock-jump conditions become 


[K+Hl 

= !W1®» 

(5.13) 

i 


(5.14) 


where the subscripts n and t refer to the normal and tangential directions, respectively. 
The tangential component of the shock-jump condition can be further simplified by use 
of Equation 5.10 to show that the jump in tangential velocity across the shock must be 
zero, i.e., 

[[Vi]] = 0 (5.15) 

Linearized Shock Jump Conditions 

The shock and wake jump conditions may also be linearized. As in the linearization of 
the field equations, a perturbation series is assumed to represent the total flow field. 
As with the moving airfoil, the motion of the shock is described by a surface vector R 
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which is composed of a mean plus an unsteady part. Hence, the surface vector as well 

as the unit vectors normal to and tangent to the shock are given by 


R(«,t) = R(s) + r(»,t) 

(5.16) 

n(s,i) = S(s) + n(*,t) = H-r^ 

(5.17) 

f (s, t) = r(s) + r (s, t) = r + 

(5.18) 

Substitution of these small perturbation assumptions for the shock displacement and 
the primitive fluid-dynamic variables into the nonlinear shock-jump conditions, and 

collection of terms of zeroth and first order gives the nonlinear mean 

« 

unsteady shock-jump conditions. The steady jump conditions are 

and the linearized 

[[Wn)\ = 0 

(5.19) 

[[rt + p\] = o 

(5.20) 

ra = o 

(5.21) 

[[(, V + H v ”]] =0 

(5.22) 

and the linearized unsteady shock conditions are 



(5.23) 

[[p + pVt+ 2pVnVn\]=0 

(5.24) 


(5.25) 

+ w.v • v + 7 1 JV.P+ l 

)*“ 


j = 0 (5.26) 


In Equations 5.22 and 5.26 the energy has been expressed in terms of the four primitive 
variables. These shock jump conditions are valid at the instantaneous position of the 
shock. This form of the shock jump conditions is not always convenient. For example, 
suppose we are trying to calculate a flow with a shock using a stationary grid. At time 
to the shock is aligned with a grid line and the shock jump conditions may be applied 
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directly. A short time later the shock will have moved a small distance away from the 
grid line. We can extrapolate the variables from the fixed grid to the new shock location 
to obtain a slightly modified set of shock jump conditions. They are given by 


[[pv n + pV n -pV t ^-^p + r^{pV n )^ = 0 (5.27) 

[[p + pVl + 2 pVnVn + r • V(P + pVZ)]\ = 0 (5.28) 

[h +v -^ +r,v H] =0 (5,29) 

[[^W. + w.v.v+ ^f. P+ \w') 

+» • • V ( (-^n P + \w) V„)]] = 0 (5.30) 


The unknowns in Equation 5.30 are the magnitude of the shock displacement r, and 
the perturbation variables p, u, v, and p on either side of the shock. Known are the 
steady flow variables p, U , V, and P. As with the linearized field equations, the time 
derivative operator ^ may be replaced by jui if the flow is harmonic. 

5.2.2 Unsteady Wake Jump Conditions 

Nonlinear Wake Jump Conditions 

Although the jump conditions presented in Section 5.2.1 were derived for the case of a 
shock discontinuity, in fact they are valid for any surface at which the flow variables are 
discontinuous including blades and wakes. The blades and wakes are surfaces through 
which mass does not flow and therefore, the flow tangency condition also applies at wake 
surfaces. However, whereas the motion of the airfoils is prescribed, the wake motion 
must be determined as part of the solution. 

Because the mass flux through the wake is zero, the conservation of tangential mo- 
mentum and energy are automatically satisfied. The remaining jump condition is the 
normal momentum. An examination of the normal momentum jump condition reveals 
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that the condition reduces to the requirement that the pressure be continuous across 
the wake. Thus, the three conditions which are applied at the wake surface are 

3ft . 


(V-n) 


upper * 


dt 


(V.n) low , r =^.n 

m = o 


(5.31) 

(5.32) 

(5.33) 


Hence, pressure and the normal component of velocity are required to be continuous 
across a wake. The tangential velocity and density will generally be discontinuous. One 
final requirement on the wake position is that the wake be attached to the trailing edge 
of the airfoil. 


Linearized Wake Conditions 

The unsteady wake conditions are linearized in a fashion similar to the unsteady shock 
jump linearization. The wake surface R is described by a perturbation series as are 
the primitive variables. These assumptions are substituted into the wake equations 
(Equations 5.31-5.33). Collection of terms of equal order provides the zeroth-order 
(mean) and the first-order (linearized unsteady) wake conditions. The steady jump 
conditions are 


(V*ff) Iower = 0 

(5.34) 

(V-ff) upper = 0 

(5.35) 

m = o 

(5.36) 

The linearized unsteady wake jump conditions are given by 


(- ■ -JT 

(5.37) 


(5.38) 

M = o 

(5.39) 
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Equations 5.37-5.39 are applied at the instantaneous wake positions. The corre- 
sponding conditions that apply at the mean wake positions are 



(5.40) 


(5.41) 

[[p + r • VP]] = 0 

(5.42) 


5.2.3 Newton Iteration Solution Procedure for Flows with Shocks 
and Wakes 

In Chapter 4, the basic Newton iteration procedure used to solve two-dimensional 
steady flow problems was discussed. The extension of the method to two-dimensional 
flows with shocks is discussed briefly below. Moretti [26,27] has pioneered the use of 
shock fitting in Euler codes. Although the idea of modelling the shock as a line discon- 
tinuity is similar, the implementation of the present shock fitting technique is different. 
Moretti uses a time-marching algorithm to determine the steady shock location. In the 
present method, however, the shock position is found using a direct method. 

The basic procedure for fitting two-dimensional steady shocks is the same as for the 
one-dimensional flows. First, an initial position of the shock is assumed. The grid is 
generated such that a double grid line straddles the shock. Hence, the flow is divided 
into an upstream and downstream region. The flow in the two regions are coupled 
together by applying the shock-jump conditions at the face centers of cells lying along 
the shock. At each of the double nodes, an additional variable is introduced which is 
the distance the shock will be displaced from this initial guess. Because shock fitting 
is part of the Newton iteration procedure, the nonlinear steady shock-jump conditions 
will be linearized by assuming that the correction to the assumed shock position is small 
as is the correction to the primitive variables. Using these assumptions, the linearized 
steady shock-jump conditions are 

[[pVn + pV n -jV t ^ + T-V(pV n )J\ = -[{pV n ]] (5.43) 

[[p + pVf + 2 WnVn - + r • V(P + pV n 2 )]] = - [[rt + p]] (5.44) 
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( 5 . 45 ) 


+ +r-V ( V,)]]=-HV,| | 

[[y>V*V„ + ?V„V -v + ^V„p+ (^~P+ \W') »» 

- kV + H ' V (feV + H v -)]] “ 

-[[(^T p + H v "]] (546) 

These boundary conditions are discretized at the centers of the cell faces lying on the 
shock using second-order accurate differencing techniques. 

Assuming that the conservation equations and the other boundary conditions have 
been discretized and assembled, Gaussian elimination is used to solve for the corrections 
to the primitive variables. Also obtained as part of the solution is the correction to the 
shock position. First, the primitive variables are updated using the corrections. Then, 
using the correction in the shock position, the grid is deformed so that the double grid 
lines coincide with the new shock location. But we are not quite done. Because the 
grid nodes have moved, the primitive variables stored at the old node locations must 
be extrapolated to the new node locations. This completes one Newton iteration. The 
entire process is repeated until convergence. 

Described above is the procedure used to fit shocks using the four conservation 
equations. If, however, the steady flow is homoenergetic, then the energy equation may 
be eliminated. The procedure for fitting shocks is essentially the same as just described 
except that there will only be three shock-jump conditions. 

Although not discussed in detail here, the procedure for fitting wakes is similar to the 
procedure for fitting shocks. The nonlinear steady wake jump conditions are linearized 
for use with the Newton iteration procedure. At each step of the iteration, the double 
grid line which represents the wake is moved to coincide with the new estimate of the 
wake location. 
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5.3 Periodic and Far Field Boundary Conditions 

Periodic Boundary Conditions 

The periodic boundaries upstream and downstream of the blades requires special treat- 
ment. The flows which we will consider are harmonic in time and periodic in the cascade 
circumferential direction. The line which extends from the leading edge of the blade 
forward in the axial direction is a periodic boundary. For steady flow, the state of the 
fluid along this line will be equal to the state of the fluid along the corresponding line 
for the adjacent blade. Similarly on the wake boundary, the flow just above the wake of 
one blade will be the same as the flow just above the wake of the adjacent blade. These 

m 

are known as the periodicity conditions. They allow us to calculate the flow through 
the entire cascade by considering a single blade passage. 

For unsteady flow, the periodicity condition is slightly different. In this case the 
unsteady flow along the periodic boundary will be the same as the flow along the periodic 
boundary of the adjacent blade but phase shifted by the interblade phase angle. For the 
flutter problem, the blades are assumed to vibrate in such a way that a travelling wave 
is set up in the rotor. The interblade phase angle a is the phase difference between the 
motion of a given blade and its neighbor. Mathematically, this periodicity condition is 
expressed as 

/(x, y + s) = /(*, y)e ja (5.47) 

where / is any one of the perturbation variables and s the distance between blades. This 
boundary condition is applied at all the boundary nodes along the upstream periodic 
boundaries. 

Far-Field Boundary Conditions 

To analyze the upstream far-field behavior of the unsteady flow requires a bit of work. 
Because we will be solving for the unsteady flow in a cascade numerically, the computa- 
tional domain cannot extend indefinitely in the axial direction. Usually, the computa- 
tional domain will extend only about one chord length in front of the blade row. At this 
upstream boundary of the computational domain, we need to apply an appropriate set 
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of conditions which will allow unsteady disturbances to pass out of the domain without 
being reflected. The analysis below is somewhat similar to the characteristic boundary 
analysis used in time marching-schemes. In such schemes (unlike the present method), 
it is often assumed that any unsteady disturbances which impinge on the upstream 
boundary strike it with the wave fronts parallel to the boundary. That is, the partial 
derivative with respect to the tangential direction is negligible. This assumption results 
is the well known one-dimensional characteristic boundary conditions. The character- 
istic variables represent vorticity, entropy, and two isentropic pressure waves. However, 
if the wave fronts are not parallel to the boundary, the wave will be partially reflected 
back into the computational domain. For steady state calculations, this slows conver- 
gence. For unsteady flow calculations, this produces errors in the predicted unsteady 
flow quantities. Recently, several investigators have suggested higher-order nonreflect- 
ing boundary conditions which produce smaller reflections [38,39,40]. Fortunately, in 
the present analysis, exact nonreflecting boundary conditions can be imposed along the 
upstream boundary. This procedure is described below. 

Consider the cascade passage shown in Figure 5.1. Note that the x- and y-axes are 
aligned with the axial and tangential flow directions. We assume that sufficiently far 
upstream of the rotor the flow is uniform. Hence the linearized Euler equations can be 
expressed as 

b -£ +b ’£ +b ’£=° (548) 

where Bi, B 2 , and B 3 are the matrices which appear in Equation 4.24, and u is the 

vector of perturbation variables [p, u, v, p]^\ If the cascade is vibrating with interblade 

phase angle <r t then the solution has a spatial period in the tangential direction of 2 xs/a. 

Because we are interested in the behavior of waves in the far-field, a more natural 

representation of the solution upstream of the rotor is given by the Fourier series 

u(x,y,t) = £ (5.49) 

• =— 00 

where Pi = (<r + 2xi)/s, and I fc< is a spatial wave number to be determined. Substitution 
of Equation 5.49 into Equation 5.48 and collection of terms leads to 

f 2 [wBi + + PiB 3 ] fi . e (M+}P.ii+iki*) = o (5.50) 

t=-oo 
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For this equation to hold, each term of the infinite series must vanish. Hence, 

[wBi + fc|B2 + u< = 0 (5.51) 

Recall that u> and fa are prescribed quantities. Equation 5.51 is then an eigenvalue 
problem for the eigenvalue Jfe< and the eigenvector u*. After some matrix manipulation, 
we can express this eigenvalue problem in a slightly more convenient form, i.e., 

- Bj 1 [u/Bi + faB s ) ifa = kifa (5.52) 

Here, the matrix on the left hand side is known, and the eigenvalues and eigenvectors 
are to be determined. These eigen-quantities are important since they determine the 
nature of the travelling waves, as well as the speed and direction in which they travel. 
In particular, it can be shown that the so-called characteristic variables are directly 
related to the eigenvectors. To demonstrate this, we make a linear transformation from 
the four primitive variables (Fourier coefficients) u, to the four characteristic variables 
w Hence we have that 


u, = Tw, (5.53) 

w, = T _1 u, (5.54) 

Substitution of Equation 5.53 into Equation 5.52 and pre-multiplication by T -1 gives 

- T _1 BJ 1 [wBi + faBs] Twj = fcjW, (5.55) 


were T is chosen so as to diagonalize the right hand side. That to say, T is the matrix of 
right eigenvectors while T _1 is the matrix of left eigenvectors. By solving this eigenvalue 
problem, it can be shown that the four characteristic variables are given by 


0 

faV + uj 

-fau 

0 

faV + OJ 

-fau 

1 

0 

0 

0 

fau 

faV + u 


yJPt(U*+V*-a?)+2PiwV+u>* 

JS 

yj0j(U->+V‘ 1 -a.‘ 1 )+10 i wV+^ 

pa 

1 

~a s 

fit 

P 


Pi 

Vi 

pi 


(5.56) 
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The corresponding wave numbers are 

U(u + PiV) + y/pfjjT 2 + V* - a 1 ) + 2p<u>V + w* 


k u = 


U 2 - a 2 


U(u + PiV) - P* {U 2 + V* - a 2 ) + 2piuV + u 2 

k2i ~ U 2 - a 2 

. ftV' + a; 

«3i = “ ■ 


* 4 ,- = - 


u 

PiV + u 

U 


(5.57) 

(5.58) 

(5.59) 

(5.60) 


The characteristic variables represent the downstream and upstream moving pres- 
sure waves, an entropy wave, and a vorticity wave for the given wave number Pi. The 
speed at which these waves propagate, and the nature of their propagation, sure related 
to their respective wave numbers. If the wave number is real, then the wave propagates 
unattenuated and is said to be superresonant. If the imaginary part of the wave number 
is negative, then the wave grows in the x-direction. If the imaginary part of the wave 
number is positive, then the wave attenuates in the x-direction. These are so-called 
subresonant modes [15]. 

Note also that the wave numbers are nonlinear functions of the temporal frequency 
u), which is an indication that the system is dispersive. This means that the energy 
stored in the waves of different wave numbers will in general travel at different speeds 
through the fluid. In a dispersive system, the phase velocity and group velocity are 
different. The phase velocity is the speed at which a wave appears to move if one 
follows the wave peak. The group velocity is the velocity a packet of waves of a given 
wave number moves. The group velocity determines the direction of movement of the 
characteristic variables. 

The phase velocity, denoted by c, is the speed at which a wave appears to move if 
one follows a wave crest. Hence, the phase velocity is given simply by 

(5.61) 

The group velocity C is the speed at which a packet of waves of a given wave number 
moves. It can be shown that this speed is given by 

C =- d J±=-( d ±Y' (5.62) 

dk \doj) v 1 


w 

c= ~k 
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If the waves propagate unattenuated, then it is the group velocity which determines 
which of the characteristic variables are to be prescribed at the far-held boundary. The 
waves which enter the domain must be specified. For the flutter problem, this means 
the incoming waves would be set to zero. For the gust response problem, the values of 
the characteristic variables would be chosen so as to model the incoming entropy and 
vorticity waves. 

Due to the square root which appears in the expression for the pressure wave num- 
bers, it is possible for the wave numbers to be complex. In this case it is the imaginary 
part of the wave number which determines which characteristic values must be specified. 
The pressure wave which decays in the x-direction is the one which must be specified at 
the upstream. Whenever the conditions are such that pressure waves do not propagate, 
the waves are said to be cut-off [41]. The condition where a wave just propagates is 
known as acoustic resonance [6]. 

Discretization of the Upstream Far-Field Boundary Condition 

The procedure for specifying the upstream boundary condition is as follows. Along the 
upstream boundary, there are J boundary points. But because the first and last points 
are related by the periodicity condition, there are really only J — 1 independent points. 
To avoid confusion, let these J — 1 points be indexed by n. Then the harmonic perturba- 
tion variables can be expressed along the far-field boundary by a linear transformation. 
So for example the pressure is given by 

iji-i 

u» = Y, ^i e30,Vn (563) 

where y n is the location of the nth node along the inlet boundary, u, is the amplitude 
of the »th Fourier component of the primitive variables, and Pi is the »th wave number. 
Equation 5.63 can be expressed in matrix form 

u = Eu (5.64) 

where now u is understood to mean all the perturbation primitive variables along the 
upstream far-field boundary, and u is understood to mean the Fourier components of 


120 


the primitive variables. Inverting this relationship gives the inverse transform 


u = E -1 u 


(5.65) 


The vector ii is composed of J — 1 submatrices denoted by u». The characteristic 
boundary conditions are applied to each of these subvectors. If the flow is subsonic at 
the inlet, then three boundary conditions must be specified for each u,-. These boundary 
conditions correspond to setting the values of the incoming pressure wave, the entropy 
wave, and the vorticity wave, and can be expressed as 


where 


K,= 



• 

K.Ui 

= f, 

(5.66) 

0 

PiV + u 

-PiU 

V/??(tf a +V^-S a )+*/»i"V>w a 

pa 
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“5* 

(5.67) 

0 

W 

PiV + u 

«S|i* 



and f, represnets the specified levels of the incoming waves. Usually, f,- will be zero. The 
exception is the case where * = 0 which corresponds to the fundamental wave number. 
For the flutter problem, fo is set to zero. For the gust response problem, fo is chosen to 
represent the incoming gust. 

The J - 1 conditions can be expressed in matrix form as 


Ku = f 


(5.68) 


where K is a matrix with 3 (J — 1) columns and 4 (J — 1) rows. This matrix is sparse 
and could be called block diagonal, with K< on the diagonals. The vector f has mostly 
zero entries except for the three entries which specify the incoming gusts. Making use 
of the inverse transform (Equation 5.65) gives 

KE -1 u = f (569) 


Note that these boundary conditions couple all of the nodes along the far-field boundary 
since KE -1 is not sparse. These boundary conditions are added to the discretized 
conservation equations by the squaring technique introduced in Chapter 4. 
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Downstream Far-Field Boundary Conditions 

The boundary conditions at the downstream far-held boundary are handled in essen- 
tially the same fashion as the upstream far-held boundary condition except for two 
important differences. First, if the outflow is subsonic, only one boundary condition 
for each tangential wave number /?,• is applied. This boundary condition corresponds to 
the upstream moving pressure wave. Secondly, the velocity (and the density) are not 
continuous across the wake. Hence, the how is divided into two parts: a continuous 
part which is continuous through the wake, and a discontinuous part which satishes 
the jump in velocity in the wake but has no unsteady pressure associated with it. The 
nonrehecting boundary condition is then applied to the continuous part only. For a 
more complete description of this technique, see Reference [37]. 

5.4 Surface Pressures 

Once the unsteady flow through a cascade has been calculated, the results can be 
used to calculate the unsteady lift and moment acting on the blade by integrating the 
surface pressure over the surface of the airfoil. This process is slightly complicated by 
the fact that the airfoil may be vibrating so that the unsteady pressures at the mean 
blade location are not the same as at the instantaneous blade location, and because 
there may be shocks oscillating on the surface of the airfoil giving rise to an anharmonic 
pressure distribution 

5.4.1 Moving Airfoils 

To determine the unsteady pressure at the surface of the airfoil requires that the 
unsteady pressure perturbation p, and the gradient of the steady pressure VP at the 
mean blade location be known. To obtain the surface pressure at the instantaneous 
blade location, one simply extrapolates from the mean to the instantaneous location. 
Hence, to first order, 

peH = peH + VP •re*'* (5.70) 

■surface I mean v 

where r is the prescribed unsteady blade motion. 
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5.4.2 The Shock Impulse 


A slightly more complicated problem is that of shock motion on the surface of an 
airfoil. Suppose now that the airfoil is stationary but a shock is oscillating on the 
surface of the airfoil. Away from the shock, the unsteady pressure is simply given 
by the perturbation pressure pe’ wt . But in regions close to the mean shock location, 
the unsteady pressure is anharmonic due to the passing of the shock. As the shock 
passes over a given point on the surface, there is an order unity jump in the pressure. 
Figure 5.3 shows this situation graphically. Plotted is the pressure on the surface 
of an airfoil at a given instant in time along with the mean pressure distribution. The 
difference in the instantaneous pressure and the mean pressure is the unsteady pressure. 
Away from the shock, the difference is first order. In the region between the mean and 
instantaneous shock locations, the pressure difference is order unity, but the difference 
is shock locations is first order. Hence, the anharmonic pressure produces first-order 
contributions when integrated to obtain the lift or moment. Because the shock motion is 
harmonic, this contribution is also harmonic to first order. For aeroelastic calculations, 
one can consider the anharmonic pressure due to the shock motion to be a harmonic 
impulse. This impulse is represented as 

Pirn pul-e*" « (5.71) 

where AP is the mean pressure jump across the shock and x, is the complex amplitude 
of shock motion along the blade. 


5.5 Summary 

In this chapter, the boundary conditions needed to analyze unsteady cascade flows 
have been presented. Of special interest are the treatments of flow discontinuities and 
far-field boundaries. Unsteady shocks and wakes are analyzed using shock and wake 
fitting procedure rather than by capturing. Fitting allows an accurate representation of 
these discontinuities without excessive grid resolution. The far-field boundary conditions 
are specified using nonreflecting boundary conditions. These boundary conditions also 
allow one to specify incoming gusts for the gust response problem. 
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Figure 5.3: Anharmonic unsteady pressure due to shock motion. Upper plot shows the 
mean and instantaneous pressures on an airfoil surface. Lower plot shows the difference 
between the instantaneous and mean pressures. Note the “impulse” due to the shock 
motion. 
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Chapter 6 


Results 


In this chapter, a variety of steady and unsteady flow solutions are presented to 
demonstrate the capabilities of the present Euler analysis. In Section 6.1, steady and 
unsteady, subsonic and transonic flows through a hyperbolic channel are examined. The 
main purpose of these examples is to demonstrate the accuracy and convergence prop- 
erties of the present method for steady flow calculations, to demonstrate the steady 
shock fitting capabilities, and to demonstrate the ability of the present method to accu- 
rately predict unsteady shock motion. In Section 6.2, the results of the present method 
are compared to a time-marching method. It is shown that the two methods are in 
good agreement so long as the level of unsteadiness is moderate, and that the present 
method is computationally much more efficient. Finally, in Section 6.3, some steady and 
unsteady cascade flows are presented. The unsteady flows are due to both blade mo- 
tion (the flutter problem) and incoming vortical and entropic waves (the forced response 
problem). Here again, the method is shown to give good agreement with semi-analytical 
and time-marching methods, but are obtained with a fraction of the computational effort 
required by the time-marching methods. 

6.1 The Flow in a Hyperbolic Channel 

The first test to be considered is the flow through a hyperbolic channel. This ge- 
ometry was studied by Emmons [42,43] in the mid 1940s and was chosen here as an 
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Figure 6.1: Typical grid used to calculate steady and unsteady flows through hyperbolic 
channel. 

initial test case for several reasons. First, the flow through this geometry is highly two 
dimensional, i.e., there is significant variation in flow properties across the width of the 
channel so that it represents more than just a slight extension of one-dimensional flow. 
Second, if a shock occurs, it will be curved because of the highly two-dimensional flow 
field. As such, the hyperbolic channel is a good test geometry on which to evaluate the 
steady shock-fitting algorithm. Finally, the published results of Emmons’ allow some 
qualitative comparisons with the present steady flow results. 

The channel itself is defined by a transformation from £,77 coordinates to x,y coor- 
dinates by 

x = sinh(£) cos(r?) (6.1) 

y = cosh(£) sin (77) (6.2) 

In the coordinate system, the channel is rectangular extending from f = —1.05 to 
f = 1.05 and 77 = 0 to 77 — 0.6. A typical computational grid used to compute the 
steady and unsteady flows through this channel is shown in Figure 6.1. The upper 
wall is curved while the lower wall is straight. Emmons calculated both subsonic and 
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transonic flows through this channel geometry using a numerical method in which a 
stream function was determined using a relaxation technique. 

6.1.1 Accuracy of the Steady Solver 

The first use of the hyperbolic flow geometry will be to check the order of accuracy 
of the steady flow solver. In Chapter 4 the claim was made that the present scheme 
is second-order accurate. In other words, if the linear dimension of the conservation 
cells is 0(Ax), then the norm of the error of the computed solution (a measure of the 
difference between the exact solution and the numerical solution) is 0(Az 2 ). Although 
the exact solution is not known for this geometry, it is known that for subsonic flows 
the total pressure should be constant throughout the channel. Hence, the norm of the 
total pressure error is also a measure of the solution accuracy. 

The following procedure was used. Grids with increasingly fine resolutions were used 
to calculate the steady flow through the hyperbolic channel. The boundary conditions 
are such that the flow is subsonic throughout the channel. Upstream the total pressure 
and density were specified to be unity. Downstream, the exit pressure was specified to 
be 0.92. The steady flow was calculated using the Newton iteration procedure outlined 
in earlier chapters. Once a converged solution was obtained, the Li norm (the root 
mean square) of the total pressure loss was calculated. 

Figures 6. 2-6. 5 show the grids used to calculate the steady flow along with the 

steady Mach number contours and total pressure contours. The grids used contained 
8 x 2, 16 x 4, 32 x 8, and 64 X 16 cells respectively. As expected, the 8x2 grid gives only 
fair results with a total pressure loss in the vicinity of the throat of about 3 percent. 
The solution is dramatically better for the 16 x 4 grid with about a 1 percent total 
pressure loss at the throat. The 32 X 8 and the 64 x 16 grid give nearly exact solutions 
with less than 0.3 and 0.1 percent total pressure loss, respectively. The Mach number 
contours for these last two cases are indistinguishable from one another. 

Shown in Figure 6.6 is the norm of the total pressure loss as a function of grid 
resolution. When plotted on a logarithmic plot, the slope of the curve gives the order of 
accuracy of the method. Note in particular that the slope of the curve is very nearly 2 
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Figure 6.2: Very coarse resolution grid used to calculate steady subsonic flow in a 
hyperbolic channel. Also shown are the computed total pressure and Mach number 
contours. 
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Figure 6.4: Medium resolution grid used to calculate steady subsonic flow in a hyperbolic 
channel. Also shown are the computed total pressure and Mach number contours. 
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Figure 6.6: Norm of total pressure error as a function of grid resolution. Note the slope 
is very nearly 2 demonstrating second-order accuracy. 

indicating the scheme is second-order accurate. 

6.1.2 Convergence of Subsonic Solutions 

In this section, the convergence properties of the Newton iteration procedure will 
be numerically investigated. Consider again the solution presented above for the 32 x 8 
grid. The converged solution is shown in Figure 6.4. To obtain this solution requires a 
finite number of iterations. In Figure 6.7, the norm of the residuals of the discretized 
Euler equations is plotted as a function of iteration number. Note the extremely fast 
convergence. This behavior is typical of the Newton iteration procedure. After four iter- 
ations, no further progress is made. Smoothing and round off errors prevent the solution 
from ever approaching the exact solution to the numerical equations. Nevertheless, this 
example demonstrates the extremely feist convergence to the final solution. 
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0,1.3 Steady Choked Flow in a Hyperbolic Channel 

In this section, transonic (choked) flow in a hyperbolic channel will be considered. 
As with the subsonic case, the total pressure and density at the inlet is specified to 
be unity. Downstream, the static pressure is set to 0.85. At this back pressure, the 
flow will be choked. The flow becomes supersonic as it passes through the throat and 
“shocks down” aft of the throat. This test case is a good one for evaluating the shock 
fitting portion of the present method. 

The shock is initially supposed to span the channel width at some somewhat ar- 
bitrary location. The grid is generated so that a double grid line is imposed at this 
shock. One iteration of the Newton procedure is run to obtain a new estimate of the 
solution along with a new estimate of the shock location. The grid is then altered so 
that the new shock grid lines are aligned with the new estimate of the shock position. 
The updated solution is extrapolated from the old grid to the new grid, and the entire 
procedure is repeated until it converges. 

Figure 6.8 shows a typical convergence history. Shown are the Mach number con- 
tours and shock positions after each iteration. Note that the shock converges rapidly to 
its final position. This is somewhat easier to see in Figure 6.9 Shown in one channel are 
the shock positions after each iteration. Note that after the third iteration, the shock 
is nearly at its final converged position. Additional iterations produce shock positions 
which virtually overplot this one. 

Figure 6.10 shows the total pressure in the channel. Note that ahead of the shock, 
where the flow is isentropic, the total pressure is constant. Downstream, however, there 
is some total pressure loss due to the shock. The contours of total density are aligned 
with the local streamlines. Note the gradient in the total pressure as one moves from 
the lower wall toward the upper wall. This gradient appears because the shock strength 
varies from one side of the channel to the other. The Mach number is higher at the top 
of the channel and hence the shock is stronger resulting in a larger total pressure loss. 
Note further that the shock is normal to the wall at both the upper and lower walls. 
Because oblique shocks turn the flow coming into the shock, the shock must be normal, 
or there must be another shock reflected off of the wall. 
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Figure 6.8: Convergence history of Mach number contours for choked flow problem. 
Note that convergence is extremely fast after the first few iterations. 


135 



divergence history of shock position. Position is converged after about 



Figure 6.10: Total pressure contours in choked hyperbolic channel. Note that ahead of 
the shock, the total pressure is constant. Downstream, lines of constant total pressure 
follow streamlines. The shock is strongest along upper channel wall. 
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6.1.4 Effect of Grid Resolution on Shock Location 


For Euler codes which use shock capturing to locate the position of the shock, 
a fine grid is needed in the vicinity of the shock to accurately position the shock. 
Recently, investigators have used adaptive gridding techniques to refine the grid around 
shocks [44,45]. Hence, fine grid resolution is used only where needed, making these 
schemes more efficient for a given level of desired accuracy. With shock fitting, however, 
the grid need not be extremely fine at the shock. To demonstrate this feature, the case 
presented in the previous section was computed on a 16 x 4, 32 x 8, and 62 x 16 grid. The 
converged computed solutions are shown in Figure 6.11. Note that even the coarsest 
grid gives good agreement with the fine grid solution. 

More interesting is the predicted shock location. Shown in Figure 6.12 are the 
computed shock positions plotted in a single channel. Note that the computed shock 
positions for the 32 x 8 grid and the 64 x 16 grid are identical to within the thickness of 
the plotted lines. Even the coarse 16 x 4 grid computed shock position agrees extren !y 
well with the other two solutions. Although the difference is perceptible, it is extremely 
small, considerably smaller than the size of the finest grid cells. This is an important 
result. Shock fitting provides an accurate estimate of the shock location even on coarse 
grids. 


6.1.5 Unsteady Subsonic Flow in a Hyperbolic Channel 

In the previous section, steady subsonic flows in a hyperbolic channel were ana- 
lyzed. These solution are needed before one can calculate the unsteady flow since the 
full nonlinear equations are linearized about the steady solution. In this section, the 
unsteady flow in the hyperbolic channel will be examined. The steady flow is that first 
presented in Section 6.1.1. The flow is everywhere subsonic. The inlet total pressure 
and density are 1.0, while the exit static pressure is 0.92. The unsteadiness is introduced 
by perturbing the upstream total pressure and density while holding the downstream 
pressure constant. For the first case, the grid is a 32 x 8 cell grid. Upstream, the 
total pressure and density variations are 1.0 and 0.714 respectively. This choice of total 
pressure variation is the linearized approximation to an isentropic total pressure and 
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Figure 6.11: Effect of grid resolution on computed transonic flow solution in hyperbolic 
channel. Note that accurate solutions are obtained with fairly coarse grids. 
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Figure 6.12: Effect of grid resolution on computed shock position. Note that the coarse 
grid gives essentially the same estimate of the shock location. 

density variation. The frequency of this variation is 1.0. 

Figure 6.13 shows the real part and imaginary part of the calculated pressure along 
the lower wall of the channel. Note that the unsteady pressure goes to zero at the exit of 
the channel as specified. Note also the amplification of the pressure at the throat of the 
channel. Unfortunately, there are no other known results with which to compare these 
results. In a later section, the present method will be used to compute the unsteady 
flow in another channel geometry. There, the present method will be compared to a 
time-marching Euler code. 

0.1.6 Unsteady Transonic Flow in a Hyperbolic Channel 

In the previous section, the unsteady subsonic flow through a channel was calculated. 
Although an interesting test case, it does not demonstrate the true strength of the 
linearized Euler method. Subsonic flows can be handled perfectly well using linearized 
full potential methods. To predict unsteady small disturbance flows with shocks, the 
linearized Euler method will give a more accurate representation of the actual flow. In 
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Figure 6.14: Shock position in channel for various levels of total inlet pressure to exit 
static pressure ratio. 

this section some unsteady transonic flows will be considered. 

Quasi-Steady Transonic Flow 

Consider the full nonlinear but steady flow through the hyperbolic channel. As the 
back pressure is decreased (or the inlet total pressure is increased), the shock will move 
further back in the channel as shown in Figure 6.14. One can then consider the shock 
position as a function of the channel pressure ratio. Figure 6.15 shows the calculated 
steady shock position along the lower wall as a function of this pressure ratio. The slope 
of this curve is then approximately the distance the shock will be displaced for a small 
change in pressure ratio. 

Returning to the unsteady theory, if the frequency of the inlet total pressure per- 
turbation is small, then the unsteady shock motion will be in phase with the inlet total 
pressure perturbation, and the amplitude of the shock motion per unit total pressure 
ratio perturbation should equal the slope of the steady shock position versus pressure 
ratio curve. This provides a consistency check of the low frequency capabilities of the 
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Figure 6.15: Shock position along lower channel wall as a function of total inlet pressure 
to exit static pressure ratio. 
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Figure 6.16: Quasi-steady pressure distribution in hyperbolic channel along lower wall 
for isentropic increase in inlet total pressure. 

unsteady code. The steady flow of Section 6.1.3 was used as the base flow for the un- 
steady calculations. Recall that the steady inlet total pressure and density were set to 
unity while the static pressure was 0.85. For the unsteady calculations, the inlet total 
pressure perturbation was set to unity while the total density perturbation was set to 
0.714. The excitation frequency was 0.0. Downstream, the static pressure was held 
constant. 

Figure 6.16 shows the real part of the “unsteady” pressure distribution in the 
channel. Of course, the imaginary part is zero. The predicted shock displacement at 
the lower and upper walls is 3.598 and 2.436 respectively. Note that the upper part of 
the shock does not move as far as the lower part of the shock. This is because the Mach 
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number varies more quickly at the upper part of the channel. Hence, the change in shock 
strength per unit distance moved is stronger along the upper wall. This has the effect 
of restricting the shock motion more along the upper part of the wall. Also, the shock 
motion is finite. As obvious as this sounds, this is an important feature of the present 
method. The motion is finite because the unsteady shock strength is correctly predicted 
(for this case the quasi-steady motion) by the linearized shock jump conditions. This 
is not true of the linearized full potential methods. Because these methods assume an 
isentropic shock, there is no mechanism to fix the location of the steady shock in choked 
flows. Therefore, as the frequency becomes small, the shock motion becomes infinite. 
Said another way, the slope of the shock position versus back pressure curve is infinite. 

As a consistency check, the slope of the shock position curve (Figure 6.15) for a 
steady inlet total pressure to exit pressure of 1.176 is 3.062. Next, it is observed that if 
the exit pressure is held constant, then 

dX ‘ = 1 dX > (6.3) 

a(P*/P«dt) Pexit dPt 

Hence, for an exit pressure of 0.85, the change in shock position per unit change in 

upstream total pressure is 3.602. The predicted shock motion using the linearized 

unsteady Euler theory at zero frequency is 3.598. This excellent agreement shows that 
the linearized unsteady theory is consistent with the nonlinear steady theory for very low 
frequencies. Next, truly unsteady (not quasi-steady) shocked flow will be considered. 

Unsteady Transonic Flow 

For this example, the unsteady boundary conditions are the same as for the quasi- 
steady example except that the excitation frequency is now 2.0. The calculated unsteady 
pressure distribution along the lower channel wall is shown in Figure 6.17. Note that at 
the exit of the channel, the unsteady pressure is zero as specified. This sort of boundary 
condition reflects incident waves. That is to say that waves moving downstream will 
reflect sending waves upstream through the channel. 

Figure 6.18 shows the unsteady shock motion. Shown is the shock position at 
four different times in one complete cycle. Note three important features of the shock 
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Figure 6.18: Unsteady shock motion at four times during harmonic cycle. Note that 
upper portion of shock leads lower part. 

motion. First, the shock remains normal to the wall at all times. This is a necessary 
condition if the flow is to remain tangent to the wall. Secondly, the upper part of the 
shock does not move as far as the lower part of the shock. The shock is stronger at 
the upper wall and must move through higher Mach number gradients. This tends to 
restrict the motion of the upper part of the shock. And third, the upper part of the 
shock leads the lower part. This is a slightly more complicated feature. Along the 
upper wall ahead of the shock, the Mach number tends to be higher than along the 
lower wall. This means that the forward moving pressure and convection waves tend 
to arrive at the shock earlier at the upper part of the shock. Similarly, downstream of 
the shock, because of the higher total pressure loss along the upper wall, the backward 
moving pressure wave moves faster along the upper wall than along the lower. Hence, 
any reflected waves will arrive at the upper part of the shock more quickly than at the 
lower part. Both of these effects tend to cause the upper part of the shock to lead the 
lower. 


6.2 Two-Dimensional Transonic Diffuser Flow 

In this section, the steady and unsteady flow in a transonic diffuser will be examined. 
Steady and unsteady flows in this geometry have been numerically modelled [46] using 
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a time-accurate time-marching code. In this section, the results of the present method 
will be compared those obtained from the time-accurate time-marching Euler code. 
The advantage of a time-marching code is that it fully accounts for nonlinear effects 
in unsteady flow. Unfortunately, the number of time-accurate time steps needed to 
calculate the unsteady flow is dominated by the smallest computational cells in the 
domain. Hence, if the grid resolution is very fine to accurately model the moving shock 
behavior using shock capturing, the computational effort required to determine the 
unsteady flow will be great. As will be shown, even moderate levels of unsteadiness can 
be predicted accurately and efficiently using the linearized Euler method. Large levels 
of unsteadiness, however, produce highly nonlinear flows which cannot be represented 
in the linearized framework and are better analyzed using a time-marching technique. 

The geometry of the diffuser is shown in Figure 6.19. Normalizing lengths by the 
throat height, the channel is 12.5 long, the inlet is 1.4114 high, and the exit is 1.5 high. 
Note in particular the gentle divergence of the channel after the throat. This will tend 
to make for large shock motions for low frequency inlet or exit disturbances. 

Two flow conditions will be analyzed using the present method and the time- 
marching flux vector splitting method of Allmaras and Giles [46], The first is simply 
a steady nominal flow through the channel. Upstream, the total pressure and density 
are specified. Downstream, the back pressure is given. For a low enough back pressure, 
a shock will form aft of the throat. The second flow to be considered is the same as 
the first except that the back pressure is varied harmonically in time. This causes the 
shock to move fore and aft in the channel. 

6.2.1 Steady Diffuser Flow 

Consider the steady flow through the diffuser. Upstream the total pressure and 
density are specified. Downstream the static pressure is specified to be 0.80 of the 
upstream total pressure. At this pressure ratio, the flow is choked and a shock forms 
downstream of the throat. The grid used to calculate the steady flow using the present 
method, and the grid used by Allmaras [47] to calculate the steady flow are shown 
in Figure 6.19. Because shock fitting is used in the present method, less resolution is 
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Figure 6.19: Grids used to calculate steady and unsteady flows in transonic diffuser: 
a) grid used with present method for both steady and unsteady flow calculations; b) 
grid used by Allmaras for steady flow calculations; and c) grid used by Allmaras for 
unsteady flow calculations. 
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needed to accurately model the shock. The grid used with the present method is 40 x 10 
cells. Allmaras used an 80 x 16 cell grid. 

The computed steady flow solutions for the two methods are shown in Figure 6.20. 
The Newton iteration procedure used in the present method converged in less than ten 
iterations. Note the generally excellent agreement between the two methods. The major 
difference between the two is the representation of the shock. Nevertheless, the fitted 
shock position found using the present method falls within the region of the smeared 
captured shock of the time-marching algorithm. Also, the total pressure loss after the 
shock agrees well for both methods. 

6.2.2 Unsteady Transonic Diffuser Flow 

In the second half of this example, unsteadiness is introduced by varying the down- 
stream pressure sinusoidally in time. The mean pressure is the same as above and the 
amplitude of the pressure perturbation is 10 percent of the mean exit pressure. The 
reduced frequency of the forced excitation is 3.125 based on the sonic speed of sound 
and the diffuser length. For the unsteady calculation made using the present method, 
the same grid was used as for the steady calculation. Allmaras used a slightly different 
grid for his steady calculation which is shown in Figure 6.19. Note that there is a region 
of increased grid resolution downstream of the throat. This area roughly corresponds to 
the region over which the shock will move. The time-marching code was marched three 
periods of excitation to achieve a periodic solution. Because the two methods were run 
on different types of computers, it is difficult to compare CPU times precisely. A rough 
estimate, however, is that the time marching code required about 50 times more CPU 
time than the linearized Euler code. 

The calculated pressure contours are shown side by side in Figure 6.21. Shown 
are the instantaneous pressure contours at eight points during a single cycle for the 
present method and the time-marching method. Note that although the contours are 
not identical, there is good qualitative agreement between the two methods. This can 
also be seen in Figure 6.22. Plotted are the instantaneous pressures along the lower 
wall at eight points in the cycle. Here it is seen again that qualitatively the two solutions 
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Figure 6.20: Steady Mach number and total pressure contours as calculated using 
present method and time-marching scheme. Note excellent agreement between the two 
solutions. 
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Figure 6.21: Instantaneous pressure contours at eight points in one cycle of unsteadiness: 
a) present method; and b) time-marching results. Note the good qualitative agreement 
if not in detail. Reduced frequency w based on sonic speed of sound and diffuser length 
is 3.125. 
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Figure 6.23: First harmonic component of unsteady diffuser flow. Note almost perfect 
agreement despite the fairly moderate level of unsteadiness in the diffuser. 

agree with one another, but in detail they differ somewhat. 

From these comparisons, one might say that the agreement is only fair. However, 
the differences between the two are dominated by the higher harmonic content which 
is predicted by the time-marching scheme, but are not accounted for in the linearized 
Euler method. The first harmonic of the time-marching solution agrees quite well with 
the linearized code. To demonstrate this, the time marching results were Fourier trans- 
formed to extract the first harmonic content. The first harmonic component of the 
solution is plotted in Figure 6.23 along with the linearized Euler solution. Notice the 
almost exact agreement between the two theories. Hence, if one is mainly interested 
in the fundamental frequency component, the linearized Euler method will give good 
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results at moderate levels of unsteadiness. 

For moderate levels of unsteadiness (say 10 percent) in the flow, the linearized 
Euler method gives an accurate description of the first harmonic content of the flow 
in the diffuser. For a higher level, however, the unsteadiness will be too large to be 
handled within the framework of linearization. As an example, consider the case were 
the variation of the back pressure is increased to 20 percent of the mean flow. The 
resulting flow as calculated by Allmaras and Giles [46] is shown in Figure 6.24. Notice 
that the shock forms toward the rear of the channel and moves forward toward the 
throat. It weakens and nearly disappears as a new shock is forming. Not only is the 
shock motion not sinusoidal, it is not even continuous. Such cases clearly cannot be 
analyzed using the linearized Euler method. 
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Figure 6.24: Diffuser flow variation due to a sinusoidal exit pressure excitation of 20% 
amplitude and reduced frequency of 3.125. The mean exit pressure to inlet total pressure 
ratio is 0.80. From [46]. 



6.3 The Gostelow Cascade 


The next geometry considered is that of a cascade of Joukowski airfoils. The 
steady incompressible flow through this two-dimensional cascade has been examined by 
Gostelow [48] using a conformal transformation technique. The unsteady incompress- 
ible flow has been analyzed by Atassi and Akai [3,4] using a semi-analytical technique. 
Therefore, since steady and unsteady solutions exist for this geometry, it is a useful 
test case with which to evaluate the present linearized Euler method. First the steady 
flow through the cascade will be calculated using the present method. Then several 
different unsteady flow problems will be examined. These include unsteady flow due to 
the torsional vibrations of airfoils, and those due to incident vortical and entropic gusts. 

6.3.1 Steady Flow Through the Gostelow Cascade 

In this section, the present method is used to calculate the steady flow through the 
Gostelow cascade. An 80 x 16 node grid was generated for a single blade passage using a 
mixed algebraic and elliptic grid generator. The initial grid for the steady calculations is 
shown in Figure 6.25. The grid generator used is based on the elliptic grid generation 
technique of Thompson, Thames, and Mastin [49]. The grid is generated using this 
technique and then, in a post processing step, the grid is modified algebraically. For 
example, the part of the grid aft of the trailing edge was algebraically generated. Note 
that the normal grid lines are straight. This simplified the wake fitting somewhat. 
Ahead of the the airfoil, the elliptically generated grid was modified by fanning out 
the streamwise grid lines algebraically so that the nodes along the upstream far-field 
boundary are uniformly spaced. This improves the performance of the upstream far- 
field nonreflecting boundary conditions. The quality of the grid will often determine the 
quality of the computed solution as much as the computational method itself. To get 
good solutions with moderate grid resolution required that grid points be clustered in 
regions of high flow gradients. Figure 6.26 shows the grid clustering around the leading 
edge. 

The boundary conditions for the steady flow problem are as follows. Upstream of 
the cascade the total pressure and density were specified to be unity. The incoming flow 
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| Figure 6.25: Initial grid used to calculate steady flow in Gostelow cascade. Note initially 

the wake is assumed to be aligned in the axial direction. 
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Figure 6.26: Expanded view of the grid around the leading edge of the Gostelow airfoil. 
High grid resolution is needed in this area to achieve good computational results. 
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angle was specified to be 53.5 degrees with respect to the axial direction, or 16 degrees 
with respect to the blade chord. Downstream, the static pressure was specified to be 
0.98. This means the downstream Mach number is 0.170. 

The Newton solver was used to calculate the steady flow. Shown in Figure 6.27 are 
the computed steady Mach number contours for this case. Note that the highest Mach 
number is about 0.34. Hence, the flow is somewhat compressible. Shown in Figure 6.28 
are the total pressure contours. Because the flow is nominally isentropic, the total 
pressure should be 1.0 throughout the domain. The maximum total pressure error 
for this case is about 0.01. Because the flow is somewhat compressible, the computed 
solution cannot be compared to the analytic solution of Gostelow [48] directly. However, 
by using the Prandtl-Glauert transformation, the computed solution can be corrected 
so that it more closely approximates the incompressible solution. The result is shown 
in Figure 6.29 along with the analytical solution. Note the good agreement between 
the two solutions. 

Because wake fitting is used in the present method, the converged grid (shown in 
Figure 6.30) is different from the initial grid. Note that the streamwise grid line starting 
at the trailing edge of the airfoil lies along the computed wake position. The computed 
outflow angle is 29.65 degrees compared to the 30.02 degrees found by Gostelow. Shown 
in Figure 6.31 is a blowup of the trailing edge region of the grid. Note that the wake 
departs the trailing edge of the airfoil at the metal angle and then turns upward slightly 
to its final outflow angle. 

6.3.2 Unsteady Flow in the Gostelow Cascade 

In this section, three types of unsteady flows will be considered. First, the unsteady 
flow induced by the motion of the airfoils will be considered. This is the flutter problem. 
The goal is to determine the unsteady pressure felt by the blades due to this prescribed 
airfoil motion. Next, two gust response problems will be considered. The first is a 
vortical gust due to an upstream obstruction. The second is an entropy gust. 
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Figure 6.28: Steady flow total pressure contours. Note total pressure is nearly constant 
everywhere except close to the airfoil and wake surface. 
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Figure 6.29: Pressure along surface of airfoil. The Prandtl-Glauert correction was used 
to convert the calculated compressible flow pressure coefficient to the incompressible 
equivalent. 
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Figure 6.30: Converged Gostelow cascade grid. As part of the wake fitting method, aft 
part of grid is modified so that trailing edge grid line is aligned with wake. Calculated 
exit flow angle is 29.65 degrees. 
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Figure 6.31: Expanded view of trailing edge region of airfoil. Note the curvature in the 
wake. 
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Aerodynamic Response Due to Torsional Vibration 

For the first example, the unsteady flow induced by a cascade of vibrating airfoils is 
considered. The airfoils are assumed to vibrate in torsion about their midchords with 
a reduced frequency u 7 of 0.40 based on chord and upstream velocity. (The reduced 
frequency is given by u; = uc/V - <».) This problem was analyzed using the present 
linearized Euler analysis. The calculated unsteady pressure distribution on the surface 
of the airfoils is shown in Figures 6.32 and 6.33. Because of the difficulty in imple- 
menting the mean flow gradient terms described in Chapter 5 needed to extrapolate the 
boundary conditions from the mean blade location to the instantaneous blade location, 
these terms were “switched off” in the present analysis. This example was also analyzed 
by Atassi and Akai [3] using a semi-analytical technique that they developed to analyze 
incompressible flows. These results are also shown Figures 6.32 and 6.33. Note the 
generally good agreement between the present method and the results of Atassi and 
Akai. 

Shown in Figure 6.34 are the real and imaginary parts of the unsteady pressure 
contours. Note, in particular, two important points. First, because of the low reduced 
frequency, the unsteady pressure is nearly real. That is to say that the flow is nearly 
quasi-steady. The real part of the flow is then primarily due to the upwash induced 
by the rotation of the airfoil. The imaginary part is due to the velocity of the airfoil 
surface which is small since it is proportional to the reduced velocity. Secondly, note 
the behavior of the real pressure contours near the upstream far-field boundary. This 
smooth behavior is due to the nonreflecting boundary conditions. 

This same case was next solved for a range of interblade phase angles. Shown in 
Figure 6.35 is the imaginary part of the unsteady moment felt by the airfoil as a function 
of interblade phase angle. If the imaginary part of the moment is positive, then work 
is done on the blade over one complete cycle. Hence, the torsional blade vibration will 
grow and flutter will occur. Note that for interblade phase angles between 0 and 100 
degrees, the present method predicts that flutter will occur. These results, although 
not in exact agreement with the results of Atassi and Akai, do show similar trends. 

Atassi and Akai later included the mean flow gradient terms in their analysis [4] and 
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Figure 6.32: Real part of unsteady pressure on the surface of a pitching airfoil. Reduced 
frequency JJj = 0.4, interblade phase angle a = 7r. Lower plot is from [3]. 
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Figure 6.33: Imaginary part of unsteady pressure on the surface of a pitching airfoil. 
Reduced frequency To = 0.4, inter blade phase angle a = n. Lower plot is from [3]. 
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demonstrated that the predicted unsteady aerodynamic loads on moving airfoils are 
strongly influenced by the large mean flow gradients at the leading and trailing edges. 
Hence, these terms should be included in the linearized Euler analysis. 

Aerodynamic Response Due to Vortical Gust 

We wish to analyze the situation shown in Figure 6.36. Upstream of the rotor there 
may be some obstructions in the non-rotating frame which cause an approximately 
sinusoidal defect in the velocity. In the non-rotating frame, the flow is parallel, the 
pressure is uniform, and is uniform. In the rotating frame, however, the flow is unsteady. 
The blades see a sinusoidally varying inlet velocity and hence, the blades will feel an 
unsteady load. The object of this analysis is to predict the unsteady pressure on the 
blades due to this inlet distortion. 

To be more specific, it is assumed that the sinusoidal variation in the axial velocity 
has a wavelength of 4 blade pitches. This corresponds to an interblade phase angle 
a of —90 degrees. Furthermore, the tangential velocity in the tangential direction in 
the non-rotating frame is zero. Hence, the tangential component of velocity V in the 
rotating frame is due entirely to the rotational speed of the rotor. The unsteadiness the 
moving blades see occurs at a frequency of u = —aV/s. For the case considered here, 
this corresponds to a reduced frequency of 1.275. Time t = 0 is defined as the time at 
which the maximum velocity deficit would arrive at the leading edge of the reference 
blade if the cascade did not alter the incoming flow. 

The unsteady pressure was calculated using the present method on the grid shown 
in Figure 6.30. The code, which was run on an Alliant FX/8 computer, required about 
8 minutes of CPU time. Shown in Figures 6.37 and 6.38 are the real and imaginary 
parts of the calculated unsteady pressure distribution. By integrating these pressures 
over the airfoil surface, one can obtain the unsteady lift and moments felt by the blade 
due to the inlet distortion. Also shown for comparison is the first harmonic component 
of the unsteady pressure as calculated by Giles using the program UNSFLO [21,22]. The 
method is based on the time-accurate time-marching Lax-Wendroff method of Ni [23] . 
Giles’ code is unique in that under the right conditions a time-space transformation 
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Figure 6.36: Steady inlet distortion in the non-rotating frame is seen as an unsteady 
distortion in the rotating frame. 
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Figure 6.37: Real part of nondimensional pressure on airfoil surface due to incoming 
vortical gust. The interblade phase angle a is — tt/2, the reduced frequency cJ is 1.257. 
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Figure 6.38: Imaginary part of nondimensional pressure on airfoil surface due to incom- 
ing vortical gust. The interblade phase angle a is — tt/ 2 , the reduced frequency oj is 
1.257. 
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allows the time accurate flow field to be represented in a single passage. For low Mach 
numbers and low interblade phase angles, however, more than one blade passage may 
be needed in the calculation. For this particular example, four blade passages were 
required. The grid was therefore fairly large with 125 X 100 nodes. Eight periods were 
calculated so as to reach a periodic state. This required about 8 hours of computer time 
on an Alliant FX/8 computer. The time history of the computed pressure distribution 
was Fourier transformed so that the results could be compared directly to the present 
method. Note that the agreement is good, although not exact. At the present time, 
it is not clear what the cause of these difference are. However, considering the good 
qualitative agreement, and the fact that the present method required a factor of 60 less 
CPU time, these preliminary results are encouraging. 

Incoming Entropy Wave 

The final example is similar to the vortical wave problem. In this case, an incoming 
entropy wave is considered. As before, the disturbance is steady in the non-rotating 
frame but unsteady in the rotating frame. Upstream in the non-rotating frame, the 
pressure and velocity are uniform. The density, however, is nonuniform with a spatial 
wavelength in the tangential direction of 4 blade pitches. Hence, in the rotating frame, 
the rotor sees a uniform velocity and pressure but an unsteady density variation. This 
sort of disturbance does not cause much of a pressure disturbance on the airfoil surface, 
but is useful for illustrating the far-field boundary conditions as well as the slip plane 
representation of the wake. 

Shown in Figure 6.39 is the real part of the unsteady density. Note that at the 
upstream fair-field boundary, the contours are aligned with the axial direction. This 
stems from the fact that the entropy wave is steady in the non-rotating frame. As 
the entropy wave is convected through the blade row, the portion which passes over 
the suction surface arrives at the trailing edge sooner that that which passes over the 
pressure surface due to the circulation around the airfoil. Hence, at the wake, the 
density is not continuous. Note that this jump in density is easily handled with the 
present wake fitting technique. Finally, the wave passes out the downstream far-field 
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Figure 6.39: Unsteady density contours for incoming entropy wave. The interblade 
phase angle a is -tt/ 2, the reduced frequency u> is 1.257. Note the jump in density 
across the wake. 
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boundary without causing any artificial reflections. 
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Chapter 7 


Conclusions and 
Recommendations 


Conclusions 

In this report, a method for predicting unsteady transonic flows in turbomachinery has 
been presented. Such an analysis is needed to model accurately the aeroelastic behavior 
of turbomachinery blading. The method is based on the linearized Euler equations 
equations and fully accounts for blade loading, blade geometry, shock motion, and wake 
motion. Unlike the linearized potential methods, the generation of vorticity and entropy 
at shocks is correctly modelled. As shown, the steady and unsteady entropy generated 
at shocks has a profound influence on shock motion, and hence, the present method is 
thought to provide significantly better predictions of unsteady transonic flows. 

The analysis is divided into two main parts. In the first part, the mean flow through 
the cascade is determined by solving the nonlinear steady Euler equations. Then, the 
nonlinear time-dependent Euler equations are linearized about this operating point to 
produce a set of linear variable-coefficient equations which describe the small-amplitude 
unsteady flow. The further assumption is made that the unsteady flow is harmonic in 
time so that the explicit time dependency is removed from the linearized equations. 
These linear equations are discretized on a computational grid and solved directly. This 
approach is significantly more efficient than the alternative time-marching technique 
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with the present method needing one to two orders of magnitude less computing time 
to solve a typical unsteady flow problem. Furthermore, because the present method 
is direct, there is no CFL restriction on the size of computational cells. Hence, grid 
resolution may be increased in local areas such as leading and trailing edges without 
much increased computational effort. 

To determine the steady flow, a Newton-iteration procedure is used which reduces 
the problem of solving the nonlinear Euler equations to one of solving a series of lin- 
earized equations, each similar to the linearized unsteady equations. The steady lin- 
earized Euler equations are discretized on a computational grid and solved directly. 
The Newton-iteration procedure, as demonstrated in numerical examples, converges 
very quickly with typical subsonic flow problems converging in about five iterations 
while transonic problems may require up to ten iterations. In addition, the similarity of 
the Newton-iteration equations to the linearized unsteady equations greatly simplifies 
the effort required to develop both a steady and an unsteady Euler code. In fact, many 
parts of the two codes are nearly identical. 

Another important feature of the present analysis is the use of steady and unsteady 
shock and wake fitting. Shocks and wakes are modelled as discontinuities with jump 
conditions applied at the shock and wake surfaces. For steady flow calculations, shock 
fitting provides an accurate prediction of the shock position without excessively fine 
grids. As demonstrated, fairly coarse grids produce surprisingly good predictions of 
shock positions. But more importantly, the accuracy of the steady solution will affect 
the accuracy of the unsteady solution. Hence, the precise modelling of shocks obtained 
using shock fitting is essential if one is to accurately predict the unsteady shock motion, 
and in turn, the unsteady aerodynamic loads on transonic airfoils. 

Results of the present method are presented for a wide variety of steady and unsteady 
flow problems. These include quasi-one-dimensional channel flows, two-dimensional 
channel flows, and two-dimensional cascade flows. In many cases, these results were 
compared to those obtained from time-marching Euler codes. Generally, the agreement 
between the two methods is quite good, except for flows with large levels of unsteadiness. 
This is especially encouraging considering the computational efficiency of the present 
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method. 

Also investigated were the effects nonlinearities have on unsteady flows. From both 
asymptotic analysis and numerical experiments, the effects of nonlinearities are seen to 
appear with the square of the flow unsteadiness. As a rough rule of thumb, the linearized 
theory correctly predicts the component of the flow at the fundamental frequency so 
long as the unsteadiness in the flow is no more than about ten percent of the steady 
flow. Hence, the linearized Euler analysis is not only useful for predicting the onset of 
flutter, but also for analyzing disturbances of moderate size such as inlet distortion and 
rotor/stator interaction. 

Recommendations for Future Work 

Although the present work has demonstrated the advantages the linearized Euler method 
offers in analyzing unsteady transonic flows, some issues need to be resolved before the 
method can be successfully applied in aeroelastic analyses of transonic cascades. The 
major issue is the fitting of shocks. The shock-fitting algorithm presented in this re- 
port requires that the shock be aligned with a computational grid line on a logically 
rectangular grid. This restricts the shock geometries which can be analyzed to normal 
shocks on fairly unskewed grids. Therefore, transonic flows through staggered cascade 
cannot be currently analyzed, with stagger. Using triangular conservation cells rather 
than quadrilateral cells may alleviate this problem. Furthermore, a shock ending in the 
flow, as in the case of a supersonic patch, cannot be modelled. Some very preliminary 
efforts have been made to solve this problem with encouraging results. 

Another more difficult problem is to include the effects of viscosity. In real flows, 
the steady and unsteady flow will be affected by such viscous effects as viscous bound- 
ary layers, shock/boundary-layer interaction, and flow separation. Attempts should be 
made to include these effects into the current inviscid models to improve their prediction 
capabilities. 
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